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INTRODUCTION 


The  feasibility  of  television  bandwidth  reduction  and  spread  spectrum  transmission  has  been 
investigated  for  a  remotely  piloted  vehicle  (RPV)  [If.  Two-dimensional  transform  methods  are 
capable  of  providing  sample  rate  reductions  by  zonal  encoding  or  bit  rate  reductions  by  combin¬ 
ing  zonal  encoding  and  block  quantization.  Differential  pulse-code  modulation  (DPCM)  has  been 
established  as  equivalent  to  transform  encoding  when  a  lower  triangular  matrix  is  used  to  trans¬ 
form  the  data  to  a  set  of  uncorrelated  signals  [2] .  Mixed  transform  coding  was  investigated  and 
found  to  be  intermediate  in  performance  between  the  two-dimensional  forms  of  transforms  used 
to  encode  the  horizontal  and  vertical  directions  of  the  image,  with  the  singular  exception  of  the 
mixed  transform  composed  of  an  orthogonal  transformation  and  DPCM.  For  this  exception, 
performance  as  good  as  optimum  (Karhunen-Loeve  transformation  and  block  quantization)  was 
achieved. 

Because  ot  the  results  summarized  above,  emphasis  was  placed  on  investigating  transforma¬ 
tions  which  could  be  easily  implemented  in  real  time  using  compact  serial  access  hardware  and 
would  achieve  as  good  a  performance  as  possible  while  minimizing  the  amount  of  memory  need¬ 
ed  at  the  transmi  ter  in  the  RPV.  A  direct  two-dimensional  implementation  of  *he  Fourier  trans¬ 
form  was  conceived,  investigated,  and  found  to  be  feasible  for  nearly  square  arrays  as  reported  in 
appendix  A.  Further  investigation  was  made  of  the  mixed  transforms  using  generalized  DCPM. 
Their  performance  in  a  noisy  channel  was  evaluated  and  found  to  be  superior  to  conventional 
DPCM  systems  as  reported  in  appendix  B.  Surface  acoustic  wave  (SAW)  and  charge  transfer 
device  (CTD)  hardware  implementations  of  the  Fou.  ier  transform  were  constructed  and  evaluat¬ 
ed  and  are  presented  in  Appendices  C  and  D.  Continued  investigation  was  made  of  the  transform 
extrapolation  techniques  and  is  reported  in  Appendix  E.  Further  development  of  the  Signal 
Processing  Interpreter  is  reported  in  appendix  F.  A  feasibility  study  of  a  row-to-column  convert¬ 
er  is  given  in  appendix  G. 


SYSTEM  ANALYSIS 


There  are  two  basically  different  approaches  to  the  redundancy  reduction  and  spread  spec¬ 
trum  transmission  of  video  images  for  a  remotely  piloted  vehicle :  ( 1 )  sample  reduction  through 
transform  encoding  and  analog  transmission;  (2)  bit  rate  reduction  through  transform  encoding 
and  block  quantization  or  differential  pulse  code  modulation.  The  former  is  the  simplest  con¬ 
ceptually,  and  the  easiest  to  implement  but  imposes  restraints  on  the  spread  spectrum  encoding 
and  on  the  transmitter.  The  latter  is  more  difficult  to  implement  but  imposes  few  restrictions  on 
the  spread  spectrum  encoding  or  on  the  transmitter  and  requires  buffer  memory  and  quantizers. 

In  either  lmplementa.-on,  transform  source  encoding  is  used  to  redistribute  the  variance 
associated  with  the  image  samples  and  then  channel  encoding  is  used  to  provide  noise  protection 
and  jam  resistance.  It  is  at  the  receiver  that  the  major  difference  between  the  two  techniques 
appears.  The  first  is  essentially  a  pulse  compression  technique,  which  is  possible  since  the  redun¬ 
dancy  of  the  image  is  removed  and  spread  spectrum  redundancy  introduced  in  its  place.  Thus, 
the  receiver  can  consist  of  a  simple  matched  filter  matched  to  the  RPV  signal,  and  no  synchro¬ 
nization  is  required  between  the  RPV  and  the  receiving  station.  The  second  implementation  is 
essentially  a  digital  communication  system  in  which  channel  synchronism  must  be  established  in 
order  to  decode  the  spread  spectrum  channel  encoding. 

The  degree  of  A  J  protection  that  can  be  obtained  using  the  digital  encoding  system  in  the 
RPV  is  not  limited  by  pulse  compression  design  considerations.  The  ability  to  opeiate  with  simple 
hard  limiting  modems  further  enhances  this  alternative.  However,  some  encoding  memory  will  be 
required  in  the  vehicle  and  synchronism  of  the  received  transmission  will  be  required  prior  to  code 
demodulation  and  image  reconstruction. 
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TRANSFORM  ENCODING 

The  capability  of  transform  encoding  to  provide  either  a  sample  rate  reduction  through 
zonal  encoding  or  a  bit  rate  reduction  through  block  quantization  is  well  known  and  document¬ 
ed.  [3]  The  application  of  these  techniques  to  RPV  image  bandwidth  reduction  and  spread 
spectrum  encoding  is  complicated  by  the  large  amounts  of  signal  processing  required  and  the 
need  for  intermediate  memory  for  two-dimensional  transforms.  In  special  cases,  these  limitations 
can  be  circumvented  and  real  time  processing  h  possible  with  fixed  impulse  response  filters. 

Previous  results  established  that  for  a  one-dimensional  data  vector  of  size  N,  the  Haar  and  the 
Fourier  transforms  could  be  implemented  v  ith,  respectively,  log2  N  +  1  and  4  real  filters. 

A  new  result  is  obtained.  The  two-dimensional  Fourier  transform  may  be  implemented  with 
four  real  filters  and  no  intermediate  memory  when  the  size  of  the  two  dimensional  transform,  N, 
is  the  product  ot  the  two  one-dimensional  transform  sizes  N  j  and  N2,  which  are  relatively  prime 
(i.e.,  N!  and  N2  have  no  common  factors).  The  size  of  the  filters  required  for  implementation  are 
then  determined  by  the  transform  size  N  and  is  either  N  if  recirculating  transversal  filters  are 
used  or  2N-1  if  recirculation  is  not  used. 

Mixed  transforms  offer  the  possibility  of  improving  the  performance  of  two-dimensional 
transform  encoding  by  combining  .he  implementation  advantages  of  a  serial  transform  in  the 
horizontal  directior  with  the  reduction  of  intermediate  memory  possible  by  the  use  of  another 
transform  in  the  vertical  direction.  Since  the  performance  is  intermediate  between  the  perform¬ 
ance  achieved  by  corresponding  2-D  transforms  of  the  separate  types,  the  Cosine  transform  is 
recommended  as  the  horizontal  transform.  The  vertical  transform  may  then  be  either  the  mini¬ 
mum  memory  implementation  of  the  Haar  transform  or  the  generalized  DPCM  Sample  reduc¬ 
tion  experiments  and  bit  rate  reduction  comparisons  indicate  that  the  combination  of  the  Cosine 
and  the  generalized  DPCM  should  be  superior  in  performance  to  other  mixed  transforms.  For 
medium  error  rate  channels,  the  combination  of  spread  spectrum  techniques  with  digital  trans¬ 
mission  should  provide  the  required  output  error  rater;. 

TWO-DIMENSIONAL  DISCRETE  FOURIER  TRANSFORM 

A  class  of  two-dimensional  discrete  Fourier  transforms  (DFT),  described  in  Appendix  A  are 
shown  to  be  equivalent  to  a  one-dimensional  discrete  Fourier  transform  whose  size  is  the  product 
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Ol  the  dimensions  of  the  two-dimensional  transforms.  I  I, is  isomorphism  allows  the  two- 
dimensional  DFT  to  be  computed  in  real  time  using  the  serial  access  implementation  of  the  Chirp 
Z  transform  (CZT)  with  four  real  filters  and  parallel  implementation  of  the  complex  arithmetic. 
No  memory  is  required  to  store  intermediate  results.  This  may  be  contrasted  with  the  case  when 
the  two-dimensional  DFT  is  implemented  by  concatination  of  one-dimensional  DFTs  and  inter¬ 
mediate  one-dimensional  transforms  must  be  stored  prior  to  transformation  in  the  orthogonal 
direction.  Unfortunately,  the  input  scan  required  for  this  two-dimensional  DFT  is  not  lexical  and 
thus  not  directly  compatible  with  conventional  raster  scanned  imagers. 

It  is  shown  in  Appendix  A  that  the  input  scan  for  the  two^iimensional  DFT  is  provided  by 
iterated  translations  on  the  torus  which  results  by  folding  the  two-dimensional  image  left  to  right 
«nd  top  to  bottom.  This  implies  that  a  linearly  addressable  imager  such  as  a  vidicon  can  be 
scanned  by  applying  appropriate  “saw-tooth"  sweeps  to  both  the  horizontal  and  vertical  deflec¬ 
tion  circuits  simultaneously.  However,  for  a  continuous  sweep  sensor  such  as  a  vidicon,  the 
continuous  output  signal  is  not  directly  the  values  of  the  two-dimensional  image,  and  the  imager 
output  must  be  sampled  at  the  appropriate  times. 

The  sampling  interval  is  determined  by  the  dimensions  of  the  two  oneOimensional  arrays 
and  can  be  interpreted  as  the  time  required  to  go  to  the  subsequent  image  point  in  the  two- 
dimensional  space  of  the  input,  which  may  not  be  adjacent  to  the  preceding  point.  For  the 
special  case  when  the  transform  sizes  differ  by  one,  then  the  subsequent  sample  point  is 
adjacent  to  the  preceding  sample  point  until  an  edge  is  encountered,  and  the  scan  is  a  raster  on 
the  torus  and  is  inclined  at  an  angle  of  approximately  45“  to  the  horizontal  axis  of  the  original 
image.  The  sweep  generation  to  achieve  the  required  sweeps  for  the  twouiimensional  array  is 
electronically  simple  and  decomposes  into  two  separate  sweep  generators,  one  for  horizontal 
sweep  and  another  for  vertical  sweep.  Each  sweep  generator  can  easily  be  implemented  as  a 

recursive  loop  where  the  initial  condition  and  sweep  increment  are  determined  by  the  particular 
blo  ;k  size  and  aspect  ratio  chosen. 

l  ora  discrete  image  sensor  such  as  a  charge  coupled  device  (CCD)  or  charge  injection  de¬ 
vice  (CID),  special  attention  must  be  given  to  the  format  of  output  from  the  image  sensor.  For 
the  CCD  array  which  has  fixed  lexical  output  scan,  an  additional  memory  or  scan  converter  is 
required.  For  a  given  two-dimensional  block  size,  this  requires  approximately  the  same  amount 
of  memory  as  is  needed  in  the  row-to-column  scan  converter  used  for  the  calculation  of  the  DFT 
by  concatination  of  the  two  one-dimensional  transforms.  For  the  CID  image  sensor,  this 


limitation  does  not  apply  and  an  independent  horizontal  and  vertical  addressing  of  the  imager  is 
possible.  Thus  for  the  CID  imager,  direct  interfacing  of  the  scan  generators  to  the  imager  is  possi¬ 
ble  and  would  allow  direct  implementation  of  the  two-dimensional  transform  with  only  three 
major  components:  ( 1 )  a  random  access  CID  array,  (2)  digital  horizontal  and  vertical  scan  con¬ 
verters,  (3)  a  serial  access  CZT  transform  of  size  determined  by  the  number  cf  points  in  the  two- 
dimensional  array. 

This  approach  is  attractive  for  small  arrays  where  the  number  of  states  required  to  address 
the  imager  is  small  and  the  size  of  the  one-dimensional  transform  is  not  too  large.  Current  tech¬ 
nology  would  probably  limit  this  approach  to  arrays  of  the  size  99  X  100  elements  or  100  X  101 
elements,  cor  large.'  arrays,  the  size  of  the  one-dimensional  transform  would  exceed  104  and  its 
implementation  by  transversal  filters  would  be  difficult.  However,  the  use  of  direct  implementa¬ 
tion  of  two-dimensional  Fourier  transforms  within  blocks  can  be  applied  to  larger  arrays  in  con¬ 
junction  with  other  processing  techniques  between  blocks. 

MIXED  TRANSFORMS 

Mixed  transforms  can  be  used  either  for  intraframe  encoding  or  interframe  encoding  de¬ 
pending  on  the  available  memory  and  the  type  of  mixed  transform  implemented.  Three  mixed 
transforms  have  been  considered  for  intraframe  a  oding:  (1)  the  discrete  cosine  transform 
(DCT)  with  the  Haar  transform,  (2)  the  DCT  with  generalized  DCPM,  (3)  small  block  size  two- 
dimensional  Fourier  with  generalized  DCPM.  In  order  for  a  mixed  transform  to  be  competitive 
with  one  of  the  conventional  two-dimensional  transforms,  it  must  offer  either  superior  perform¬ 
ance  or  simplicity  of  implementation. 

Symmetrized  transforms  which  have  cosine  basis  vectors  have  this  property.  Of  these  trans¬ 
forms,  the  odd  length  cosine  transform  appears  superior  in  implementation  since  it  can  be  imple¬ 
mented  as  the  re-;!  part  of  a  CZT  and  since  the  transform  samples  are  real  and  are  the  samples  of 
an  auto  correlation  function  which  may  then  be  extrapolated  by  well  known  techniques.  In 
simulations  of  transform  performance,  the  cosine  transform  has  been  shown  to  closely  approxi¬ 
mate  the  behavior  of  the  “optimum”  Karhunen-Loeve  transform.  All  of  these  conditions  suggest 
that  the  cosine  transform  should  be  considered  as  the  horizontal  transform  in  a  mixed  transform 
implementation. 

The  minimum  memory  requirement'  of  the  Haar  transform  and  its  implementation  by  both 
transversal  filter  and  “fast”  digital  algorithms  with  a  minimum  of  computational  steps  makes  this 
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transform  attractive  as  a  vertical  transform  for  use  with  the  cosine  transform.  In  some  applica¬ 
tions,  digital  implementation  of  the  two-dimensional  Haar  transform  are  attractive.  In  analog 
systems  using  mixed  transforms  where  sample  rate  reduction  precedes  analog  transmission,  the 
zonal  filter  used  for  coefficient  retention  must  be  asymmetric  with  less  reduction  made  in  the 
Haar  coefficients  than  in  the  other  transform  coefficient. 

The  benefits  of  mixed  transformation  implementation  and  minimum  memory  may  be 
achieved  for  digital  transmission  by  combining  a  one-  or  two-dimensional  transform  with  general 
lzed  DPCM  in  a  hybrid  system.  The  basis  of  operation  of  the  hybrid  transform  is  that  the  trans¬ 
form  decorrelates  the  image  within  its  restraints  of  transform  type,  unitary  dimensionality,  and 
block  size,  while  the  generalized  DPCM  removes  the  correlation  between  transform  blocks.  This 
hybrid  system  is  particularly  attractive  for  RPV  application  since  its  performance  is  approximate 

ly  as  good  as  the  “optimum”  Karhunen-Loeve  transform  and  its  implementation  requires  mini¬ 
mum  memory. 
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HARDWARE  IMPLEMENTATION 


When  the  orthogonal  transformation  used  in  an  image  processing  application  is  implementa- 
ble  by  transversal  filter  techniques,  real  time  operation  may  be  achieved  with  small,  low-power, 
lightweight  hardware.  For  moderate  accuracy  transforms  (e.g.  7  to  10  bits  equivalent  input 
accuracy),  sampled  data  analog  transversal  filters  such  as  SAW  and  CTD  are  attractive. 

To  investigate  the  feasibility  of  these  techniques  both  SAW  and  CTD  Fourier  transforms 
were  constnicted  and  evaluated.  The  SAW  transform  was  a  32  point  discrete  Fourier  ti  nsform 
(DFT)  implemented  at  IF  with  parallel  arithmetic  on  an  ST-cut  quartz  substrate.  The  CTD  trans¬ 
form  was  a  100  point  modified  even  cosine  transform  implemented  at  video  with  parallel  arith¬ 
metic  using  bucket  brigade  devices  (BBD).  Both  implementations  were  successful.  Spectrum 
analysis  was  performed  using  both  implementafions  .  !d  is  reported  in  appendices  C  and  D. 

Comparing  the  performance  capabilities  of  the  two  implementations,  the  development  of 
CCD  transform  implementation  is  recommended  for  the  RPV  image  bandwidth  reduction.  The 
CCD  will  allow  operation  directly  at  video  frequency  with  smaller  size  than  a  comparable  SAW 
implementation;  it  also  achieves  better  mass  production  capabilities  and  lower  production  costs. 

SAW  devices  provide  an  existing  technology  by  which  real  time  video  processing  via  the  CZT 
c  m  be  implemented.  A  32  point  complex  DFT  implemented  with  parallel  arithmetic  and 
aperiodic  convolution  is  described  in  Appendix  C.  Although  filters  with  5 1 2  taps  could  be  con¬ 
structed  on  quartz  using  large  area  photolithography,  a  CZT  of  50  microseconds  of  video  would 
require  a  delay  line  in  excess  of  300  mm.  This  size  could  be  reduced  to  something  in  excess  of 
150  mm  by  substitution  of  Bi12Ge02  as  the  piezoelectric  substrate.  However,  inis  material  is  not 
isopaustic,  ard  delay  time  compensation  is  required.  Another  consideration  when  using  SAW 
devices  for  video  processing  is  that  these  devices  are  band-pass  devices  which  need  to  operate  at 
an  intermediate  frequency  (IF)  and  IF  amplifiers  and  detectors  must  be  used  with  the  SAW  de¬ 
vice  in  order  to  achieve  video  operation. 

If  the  transmission  rate  of  the  processed  v.deo  is  less  than  the  real  time  video  input  and 
minimum  memory  is  a  consideration,  then  transforms  of  smaller  size  than  256  or  5 12  should  be 
considered.  A  32  point  transform  offers  a  good  compromise  between  transform  efficiency  which 
increases  with  increasing  transform  size  and  coefficient  storage  memory  which  increases 


discontinuously  from  a  fractional  line  store  memory  to  a  frame  store  memory  when  the  coeffi¬ 
cient  computation  rate  exceeds  the  coefficient  transmission  rate.  Although  a  frame  store  mem¬ 
ory  is  possible  using  recirculating  SAW  delay  lines,  this  much  memory  can  be  organized  more 
efficiently  by  the  addition  of  specialized  control. 

With  the  availability  of  analog  frame  storage  memory,  it  is  possible  to  compute  the  two- 
dimensional  Cosine  or  Fourier  transform  by  combination.  If  the  DFT  of  the  video  is  computed 
and  stored  row  by  row  for  a  whole  frame,  then  converted  to  column  by  column  format  and  the 
DFT  computed  on  these  columns,  then  the  two-dimensional  Cosine  or  Fourier  transform  results. 
The  feasibility  of  the  required  row  to  column  transformations  using  Bi,2Ge02  recirculating  delay 
lines  was  investigated  at  Stanford  University  and  is  reported  in  Appendix  G. 

Under  some  conditions  it  is  possible  to  compute  the  two-dimensional  Fourier  transform 
directly  using  a  one-d  r.^ns-onal  CZT  of  large  size  without  the  need  of  row-to-column  conversion 
or  other  intermediate  memory.  If  the  block  size  of  the  two-dimensional  area,  N,  is  the  product 
of  two  relatively  prime  numbers,  N,  and  N-  ,  then  an  input  scan  can  be  found  which  order's  the  N 
tv.  c-cfimsrsional  data  points  so  that  they  may  be  transformed  by  a  one-dimensional  transform  of 
size  N.  This  configuration  is  described  in  Appendix  A.  Unfortunately,  the  two-dimensional 
cosine  transformation  can  not  be  computed  directly  using  this  reordering,  and  modification  of 
the  video  sensor  is  required  in  order  to  provide  the  required  input  scan.  This  technique  appears 
attractive  for  intermediate  block  sizes  in  conjunction  with  another  transform  such  as  DPCM. 
Under  these  conditions  the  transform  efficiency  of  two-dimensional  encoding  can  be  com¬ 
bined  with  moderate  one-dimensional  CZT  implementation. 

A  BBD  with  200  taps  was  used  to  construct  a  CZT  approximation  to  a  cosine  transform. 

This  component  was  developed  by  Texas  Instruments  for  Rome  Air  Development  Center  and 
consists  of  a  cosine  and  sine  BBD  on  a  common  silicon  substrate.  The  tap  weighting,  however, 
is  not  correct  for  an  exact  cosine  transformation,  and  thus  the  circuit  was  used  in  a  spectrum 
analysis  inode  as  described  in  Appendix  D.  These  devices  operate  at  a  clock  rate  of  several 
hundred  kilohertz  and,  thus,  arc  suitable  for  slow  scan  television  processing. 
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REPORT  SUMMARY 

A 

This  report  describes  (he  progress  on  the  second  phase  of  a  NUC  program  on  image  bandwidth 
reduction  for  application  to  the  ARPA  RPV  problem  of  sending  television  images  over  spread 

spectrum  channels.  This  report  details  the  first  quarter  of  a  hardware  development  and  imple- 
mentation  phase, 

Hardware  has  been  constructed  to  perform  linear  transformations  on  the  television  signai.  Both 
surface  acoustic  wave  and  charge  transfer  devices  were  built  and  tested.  It  was  decided  to  proceed 
with  charge  transfer  devices  for  the  final  design.  The  study  of  performance  tradeoff  as  a  function 
of  block  size  and  transform  has  continued.  It  has  been  determined  that  a  Cosine  transform  to¬ 
gether  with  a  DPCM  transform  in  the  second  dimension  offers  near  optimal  performance. 

A  principal  result  of  this  quarter  is  the  determination  that  the  DPCM  and  Cosine  mixed  trans¬ 
form  are  an  optimal  performance  system.  A  second  result  is  that  charge  transfer  device  hardware 
has  been  built  which  approximates  the  final  prototype  hardware. 
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REPORT  SUMMARY 


Th,s  report  describes  the  progress  on  the  second  phase  of  a  NUC  program  on  image  bandwidth 

reduction  for  application  to  the  ARPA  RPV  problem  of  sending  television  images  over  spread 

spectrum  channels.  This  report  details  the  first  quarter  of  a  hardware  development  and  imple- 
mentation  phase. 

Hardware  has  been  constructed  to  perform  linear  transfonnations  on  the  television  signal.  Both 
surface  acoustic  wave  and  charge  transfer  devices  were  built  and  tested.  It  was  decided  to  proceed 
with  charge  transfer  devices  for  the  final  design.  The  study  of  performance  tradeoff  as  a  function 
of  block  size  and  transform  has  continued.  It  has  been  determined  that  a  Cosine  transform  to¬ 
gether  with  a  DPCM  transform  in  the  second  dimension  offers  near  optimal  performance. 

A  principal  result  of  this  quarter  is  the  determination  that  the  DPCM  and  Cosine  mixed  trans¬ 
form  are  an  optimal  performance  system.  A  second  result  is  that  charge  transfer  device  hardware 
has  been  built  which  approximates  the  final  prototype  hardware. 
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A  TWO-DIMENSIONAL  DISCRETE  FOURIER  TRANSFORM  ARCHITECTURE 


by 

J.  M.  SPEISER 
H.  J.  WHITEHOUSE 

Sensor  and  Information  Technology  Department 


Naval  Undersea  Center 
San  Diego,  California  92132 


Many  problems  of  image  processing  and  other  processing  of  two-dimensional 
data  require  tiie  rapid  generation  of  two-dimaisional  Fourier  transforms  in  a 
form  suitable  to  subsequent  electronic  processing  or  utilization.  This  paper 
describes  a  two-iimensional  transform  architecture  which  is  suitable  for  real¬ 
time  image  processing  and  which  may  be  implemented  with  sity.\11,  light-weight 
special  purpose  hardware. 


INTRODUCTION 


IWo-dtensioml  Fourier  transforms  are  ordinarily  generated  either  optically 
or  through  the  use  of  a  digital  oaiputer.  Both  of  these  methods  suffer  fran 
severe  limitations  in  many  signal  processing  applications.  A  general  purpose 
digital  ocmputer  large  enough  to  perform  tw-dirensional  Fourier  transform 
in  real  time  for  many  signal  processing  applications  is  large,  heavy,  «,>ensive, 
and  consumes  a  large  amount  of  electrical  power.  Coherent  optical  two- 
dimensional  Fburier  transform  inplaientations  are  limited  by  their  input-output 
media  and  interface  poorly  with  electronic  systems. 


A  third  method  for  performing  two-dimensional  Fourier  transforms  is  to 
perform  a  row  by  row  one-dimensional  transform  with  serial  access  hardware 
such  zs  that  described  previously  [1-3],  store  the  one-dimensional  row  transforms, 
and  perform  a  oolumn  by  oolumn  one-dimensional  transform.  This  method,  however, 
require-3  311  auxiliary  memory  to  store  the  care-dimensional  transforms  of  the 
entrure  data  field.  The  auxiliary  memory  tends  to  be  expensive  in  terms  of 
cost,  weight,  and  power  consumption,  it  is  therefore  desirable  to  seek  a 
method  of  using  the  one-dimensional  serial-access  Fourier  transform  device  with 
a  nunxrmm  of  auxiliary  hardware  to  perform  ta*>-dimensKxral  Fburier  transforms. 
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A  MHTOBIONAL  DFT  ARCHITECTURE  USING  ONE-DIRENSIONAL  HARDWARE 


The  t^iAD-dirnensional  discrete  Fnirier  transform  (DPT)  system  will  use  an 
input  scanning  device  and  a  one-dimensional  discrete  Fourier  transform  devioe 
as  shown  in  Figure  1.  Th?  two-dimensional  transform  block  size  N..  by  N„  is 
chosen  such  that  and  N2  are  relatively  prime  integers  (i.e. ,  they  must  have 
no  canton  divisor) .  The  one-dimensional  Fourier  transform  device  has  a  block 
length  of  N  =  N^N_. 

The  purpose  of  the  input  scanning  device  is  to  so  order  the  input  data  that 


the  one-dimensional  Fourier  transform  of  the  length  N^N  serial  data  string  is 


identical  to  an  N  by  N  two-dimensional  Fourier  transform  of  the  N  by  N 
12  12 

input  data  samples.  If  desired,  an  output  scanning  device  may  also  be  used  to 

provide  the  transform  output  points  in  normal  order. 

The  required  scan  may  be  derived  from  Preisendorfer's  [4’  representation 

of  a  one-dimensional  discrete  Fourier  transform  matrix  as  a  direct  product 

matrix.  Preisendorf er  shows  hew  the  one-dimensional  DPT  in  equation  (1)  is 

equivalent  to  tlie  two-dimensional  DPT  in  equation  (2)  when  N  =  N  N  and  N  ,  N 

12  12 

are  relatively  prime. 


/\ 

f 

q 


N  1 


=  £ 


p=0 
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l_Nl  N2J 


jr° 


for  j^,k^  =  0,  1, 
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. .  . . 


In  order  to  make  the  tvro  transforms  equivalent,  it  is  necessary  to  find 
pair  of  one-to-one  functions  P(jr  j2)  and  Q(k1,  k2)  such  that 


PQ  _ 
N 


jlkl  +  j2k2 


(Modulo  1) 


or 


"  hkl®2  +  :>2k2Nl  (Modulo  ttjty 


This  may  be  acoomplished  by  letting 

P^l'  ^  =  -^iN2  +  -^2N1  (Modulo  N) 

Q*kl'  =  klUlN2  +  k2U2Nl  (Modulo  N) 
where  the  constants  and  U2  are  the  solutions  of 

n2ui  =  1  (Modulo  N^) 

n^U2  =  (Module  N2) 

Equations  (7)  and  (8)  will  have  solutions  if  and  only  if  ^  and  n2  are  nutually 
prime.  Under  this  condition,  the  mappings  described  by  equations  (5)  and  (6) 
will  satisfy  the  requirement  of  equation  (3). 

Equation  (5)  may  be  interpreted  as  a  prescription  for  building  the  input 
scan  generator.  As  ^  ranges  from  0  to  N^l  and  j2  ranges  from  0  to  n  -1,  P 

will  take  all  values  from  0  to  N^-l.  These  may  be  tabulated  and  rearranged 
in  increasing  order. 


i.e.  0  =  PCj^O),  j2(0) ) 

1  =  P(jT  (1) ,  j,(l)) 

l 

l 

l 

N-l  =  Pfj-^N-l),  j2(N-l)) 

Itie  values  j^(0),  .  .  . ,  j^(N-l)  and  j2(0),  •  •  •#  j2(N-l)  nay  be  stored 
in  a  pair  of  read-only  memories  to  provide  the  input  scan  generator. 

To  eliminate  the  ROM's,  we  rote  that 

P^d)  +  jx (s) ,  j2(l)  +  j2(s))  =  p^ts),  j2(s))  tp^d),  j2 (1)) 

Therefore 

j^s+l)  =  jx (s)  +  jx  (l)  (Mod  Nx) 

j2(s+l)  =  j2(s)  +  j2(l)  (Mod  N2) 

So  that  only  simple  recursive  digital  filters  are  needed  for  the  input  scan 
generator  of  Figure  2.  With  only  a  change  in  the  prestored  constants,  replacing 
j!<l)  by  ^(1)  and  j2 (1)  by  k2(l),  tne  circuit  of  Figure  2,  will  also  function 
as  the  output  scan  generator.  Representative  input  and  output  scans  are  shown 
in  Figures  3  and  4. 

A  typical  application  of  the  two-dimensional  Fourier  transform  device  would 
be  in  a  reduced  redundancy  television  system,  as  shown  in  Figure  5.  The 
selection  table  and  selection  gate  may  be  used  so  that  only  selected  spatial 
Fourier  transform  coefficients  will  be  transmil  .ted  through  the  channel. 

In  addition  to  eliminating  the  intermediate  transform  memory,  it  is  also 
possible  to  eliminate  the  ROM's  used  in  the  CZT  implementations  of  the  one¬ 
dimensional  DFT, 


" J-  ^  WW  '««  L»-  4  u  ^1, 1.J^I  1^  _L.  ,uA  IV.  1 JLJ  1^  .t»j|U  JW1JWP.JJI  y  IP.  UUL1.  iWJ  J^ii,  ,MJ  .^_1T  lii#! 


CONCLUSIONS 

As  contrasted  with  previous  digital  fco-diirensional  Fburier  transform 
devices,  the  device  debited  here  can  be  extremely  fast  since  a  large  number 
of  multipliers  are  effectively  utilized  in  the  oorrelator  or  convolver  of  the 
one-dimensional  DPT  it  uses.  As  contrasted  with  previous  optical  two- 
dimensional  Fourier  transform  devices,  the  device  described  here  interfaces 
extremely  well  with  electronic  signal  processing  systems,  and  lends  itself 
to  potential  mass  production  via  integrated  circuit  fabrication  techniques, 
representing  a  potential  very  large  savings  in  cost  for  quantity  production. 

The  one-dimensional  EFT  which  is  utilized  may  be  inplenented  using 
acoustic  surface  waves,  charge  coupled  devices,  digital  LSI  correlators,  or 
other  serial  access  transversal  filters  or  correlators. 

Transforms  of  higher  dimensionality  may  be  implemented  by  adding 
additional  sections  to  the  scan  generator.  The  only  fundamental  constraint 
is  that  the  size,  N,  of  the  one-dimensional  DPT  device  of  Figure  3  must  be  the 
product  of  the  individual  dimension  sizes  t^,  .  ,  Np,  and  evrry  pair  NR,Ng 
must  be  mutually  prince. 
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Figure  1 .  Two-dimensional  discrete  Fourier  transform  device. 
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Figure  3.  Representative  Input  Scans 
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Figure  3.  Representative  Input  Scans 
(Continued) 
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Figure  3.  Representative  Input  Scans 
(Continued) 
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Figure  4.  Representative  Output  Scans 
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Figure  4.  Representative  Output  Scans 
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hybrid  coding  of  pictorial  data* 

by 

A.  Habibi 

Department  of  Electrical  Engineering 
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Los  Angeles,  California  90007 


Two  hybrid  codjng  systems  utj)iz1ng  ,  cascade  flf  ,  unjtary  transfomatjon 

systems  are  proposed.  Both  systems  encode  the  transformed  data  by  a 
d  „f  DPCM  systems.  The  fjrst  system  usa$  a  one.dten$ionai  ^ 

a td  where  the  second  one  emp,oys  two-dimensiona,  translations.  Theoretical 
ts  for  Markov  data  and  experimental  results  for  a  typical  picture  is  pre¬ 
sented  for  Hadamard.  Fourier,  Cosine,  Slant,  and  the  Karhunen-Loeve  transforms- 

on  1SUal  e'feCtS  °f  Channel  6rr0r  3nd  also  the  of  noisy  channel 

rati  °rmanCe  °f  ^  hybnd  SySt6m’  measured  in  terms  of  the  signal  to  noise 

o  of  the  encoder,  is  examined  and  the  performance  of  this  system  is  compared 
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I.  Introduction 


In  search  for  efficient  techniques  of  transmitting  pictorial  data  over 
digital  communication  channels  has  led  various  researchers  to  a  common  approach 
to  the  problem.  Briefly,  this  approach  is  processing  the  correlated  data  (images) 
to  generate  a  set  of  uncorrelated  or  as  nearly  uncorrelated  as  possible  set  of 
signals  which  in  turn  are  quantized  using  a  memoryless  quantizer.  The  quantized 
signal  is  then  encoded  using  either  fixed  or  variable  length  code  words  and  is 
transmitted  over  a  digital  channel.  This  is  the  general  approach  taken  in  de¬ 
signing  differential  pulse  code  modulators  (DPCM)  and  the  techniques  that  use 
unitary  transformation  and  block  quantization  as  well  as  many  other  techniques 
developed  in  recent  literature  [1-6].  Both  DPCM  and  transform  coding  techniques 
have  been  used  with  some  success  in  coding  pictorial  data.  A  study  of  both 
these  systems  has  indicated  that  pach  technique  has  some  attractive  characteristics 
and  some  limitations.  The  transform  coding  systems  achieve  superior  coding  per¬ 
formance  at  lower  bit  rates,  they  distribute  the  coding  degradation  in  a  manner 
less  objectionable  to  a  human  viewer  and  show  less  sensitivity  to  data  statistics 
(picture-to-picture  variation)  and  are  less  vulnerable  to  channel  noise.  On  the 
other  hand,  DPCM  systems  when  designed  to  take  advantage  of  spatial  correlations 
of  the  data  achieve  a  better  coding  performance  at  a  higher  bit  rate,  the  equip¬ 
ment  complexity  and  the  delay  due  to  coding  operation  is  minimal,  and  the  system 
does  not  require  the  large  memory  needed  in  the  transform  coding  systems.  Per¬ 
haps  the  most  desirable  characteristic  of  this  system  is  the  ease  of  design  and 
the  speed  of  the  operation  that  has  made  the  use  of  DPCM  systems  in  real  time 
coding  of  television  signals  possible.  The  limitations  of  this  system  are  the 
sensitivity  of  the  well-designed  two-dimensional  DPCM  systems  to  picture  statis- 
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tics  and  the  propagation  of  the  channel  error  on  the  transmitted  picture. 

In  this  paper,  two  coding  systems  that  use  a  cascade  of  unitary  transfor¬ 
mations  and  a  bank  of  DPCM  systems  are  proposed.  These  systems  combine  the 
attractive  features  of  both  transform  coding  and  the  DPCM  systems  thus  achiev¬ 
ing  good  coding  capabilities  without  many  of  the  limitations  of  each  system 
The  first  system  proposed  here  exploits  the  correlation  of  the  data  in  the  hor¬ 
izontal  direction  by  taking  a  one-dimensional  transform  of  each  line  of  the  pic¬ 
ture.  then  it  operates  on  each  column  of  the  transformed  data  using  a  bank  of 
DPCM  systems.  The  DPCM  systans  quantize  the  signal  in  the  transform  domain 
Where  It  takes  advantage  of  the  vertical  correlation  of  the  transformed  data  to 
reduce  the  coding  error.  The  unitary  transformation  Involved  is  a  one-dimen- 
Siona,  transformation  of  individual  lines  of  the  pictorial  data,  thus  the  equip¬ 
ment  complexity  and  the  number  of  computational  operations  is  considerably  less 
than  what  is  involved  in  a  two-dimensional  transformation.  The  visual  effect  of 
channel  error  on  the  encoded  pictures  and  its  effect  on  the  performance  of  the 
system  measured  in  terms  of  the  signal-to-nolse  ratio  of  the  encoder  using  a 
brnary  symmetric  channel  is  examined.  The  system  shows  little  degradation  at 
small  or  moderate  levels  of  channel  noise,  however  its  performance  is  degraded 
significantly  at  high  levels  of  channel  error.  The  system  Is  particularly 
attractive  in  the  sense  that  the  principle  can  be  expanded  to  utilize  interframe 
coding  of  television  signals.  Such  a  coding  system  would  start  by  taking  a 
two-dimensional  transformation  of  each  frame  of  the  television  signal,  then  it 
would  encode  the  transformed  signal  in  the  temporal  direction  by  a  number  of 

parallel  DPCM  encoders,  thus  exploiting  the  correlation  of  data  in  temporal  as 
well  as  spatial  directions. 
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The  second  system  proposed  uses  a  two-dimensional  unitary  transformation 
on  the  pictorial  data  divided  into  small  blocks.  The  elements  of  each  block  in 
the  transformed  domain  are  ordered  in  a  one-dimensional  array  and  are  coded  by 
a  bank  of  DPCM  systems.  A  small  block  size  reduces  the  number  of  arithmetic 
operations  needed  to  obtain  the  transformed  data,  but  it  reduces  the  effeciency 
of  the  transformed  system  since  the  elements  of  various  blocks  remain  correlated 
in  the  transformed  domain.  However,  DPCM  coding  of  the  elements  in  the  trans¬ 
formed  domain  take  advantage  of  this  correlation  and  improves  the  performance 
of  the  system. 

I i •  One-dimensional  Transformation  and  DPCM  Coding 

In  the  system  proposed  here  the  pictorial  data  is  scanned  to  form  N  lines 
with  appropriate  vertical  resolution,  then  each  line  is  sampled  at  a  Nyquiflt 
rate.  The  sampled  image  is  then  divided  into  arrays  of  M  by  N  picture  elements 
u(x,y),  where  x  and  y  index  the  rows  and  the  columns  in  each  individual  array 
such  that  the  number  of  samples  in  a  line  of  image  is  an  integer  multiple  of  M. 
One-dimensional  unitary  transformation  of  the  data  and  its  inverse  are  modeled 
by  the  set  of  equations 

M 

ui(y)  *  £  ufx.y^x)  1=1,2 . M  (1) 

y=l,2 . N 

M 

u(x,y)  =  Z  u.(y)$.(x)  (2) 

where  ^ (x)  is  a  set  of  M  orthonormal  basis  vectors.  The  correlation  of  the 
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transformed  samples  u^y)  and  u^ytt)  is  given  by 


M  M 

Ci^T'  Z  R(x.x,y+r  )M*)Mx)  m 

x=l  x=l  ’  1  '  ' 

where  R(x,x,y,y)  is  the  spatial  autocovariance  of  the  data. 

Note  that  this  equation  indicates  (1)  the  correlation  of  samples  in  each 
column  of  the  transformed  array  is  directly  proportional  to  the  correlation  of 
sampled  image  in  vertical  direction,  (2)  the  correlation  of  samples  in  various 
columns  of  the  transformed  array  is  different.  Thus,  a  number  of  different 
DPCM  systems  should  be  used  to  encode  each  column  of  the  transformed  data.  The 
block  diagram  of  the  proposed  system  is  shown  in  Figure  1.  A  replica  of  the 

original  image  u*{x,y)  is  formed  by  inverse  transforming  the  coded  samples, 
i .  e . , 

n 

u*{x,y)  =  Z  v.  (y)$>.  (x) ,  n<M  (4) 

i=l  i  1  v  ' 

The  mean  square  value  of  cooing  error  Is 

2  i  N  M  „ 

e  =  E  MN  \  [u(x,y)-u*(x,y)r  (5) 

y=l  x=l 

Using  equations  (1)  to  (4)  and  assuming  that  q.(y),  the  quantization  error 

encountered  in  the  1th  DPCM  system,  is  uncorrelated  with  u^y)  the  coding 

2  .  1 
error,  e  is 
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(6) 


N  n 


e  «  R(0, 0,0,0)  -  Z  C.{0)  +  E  ^  j  z  [u, (y)-v, (y)]2 
1=1  m  y x=l  1  1 

where  the  first  two  *erms  are  introduced  because  of  using  n  (rather  than  M) 
DPCM  systems.  Study  of  DPCM  systems  has  shown  that  [1] 


E  Cu1(y)-v1(y)]2  =  E  q|(y)  =  K(m. )e2  (7) 

2 

where  e^  is  the  variance  of  the  differential  signal  in  the  ith  DPCM  and  K(m.) 

is  the  quantization  error  of  a  variate  with  a  unity  variance  in  a  quantizer 
mi 

with  (2)  1  levels. 

Analysis  of  Lloyd-Max  quantizer  has  shown  that  K(m- )  can  be  approximated 
farily  accurately  by 

-a  m. 

K(mi)  s  b  e  (8) 

where  the  best  fit  for  a  Gaussian  variate  is  obtained  using  a  =  0.5  £n  10  and 
b  =  1.0  [4,5,7].  Study  of  other  quantization  techniques  has  indicated  similar 
results  for  various  probability  density  functions  [8].  Figure  2  shows  the 
functional  form  of  K(m^ )  and  the  accuracy  of  approximation  (8)  using  a  =  0.5 
£n  10,  b  =  2.0  for  a  random  variable  with  a  double-sided  exponential  probability 
density  function  using  an  instaneous  companding  quantizer  [9]. 

From  published  results  [1],  [6]  the  variance  of  the  differential  signal  in 
a  DPCM  system  with  the  m*^  order  linear  predictor  is 


O  Ml 

ei  ’  M°>  -  f,  Vf 
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where  A(j  are  related  to  C,(j)  by  m  algebraic  equations. 
Substituting  (7)  and  (8)  In  (6),  e2  Is 


2  n  ,  n  "ami 

e  =  R(0, 0,0,0)  -  L  C  ■  (0)  +  It  Z  e  e2 

1=1  1  M  i=i  1 


(10) 


where  the  error  is  defined  in  terms  of  n  and  rr^ ,  1»i 

continuous  variables  and  minimizing  e2  with  a  constraint  E 

i=l 


Mh  i 

m.  =  —  +  1 
i  n  a 


An  e. 


1 


n 


'  n  2 

n  i=l 


An  e 


Treating  m..  as 
mi  =  W1^  give 


(11) 


where  n  is  chosen  such  that  c2  is  minimum*  and  the  quantizer  in  the  ith  DPCM 
system  will  have  (2)  1  levels,  m.  as  obtained  from  equation  (11)  is  modified 
as  discussed  in  [4]  and  [5].  Note  that  expressions  (10)  and  (11)  are  similar 
to  those  obtained  for  coding  the  transformed  data  by  memoryless  quantizers, 
the  difference  being  that  here  the  variance  of  the  differential  signal  rather 
than  the  variance  of  the  transformed  data  is  used. 


Figure  3  shows  the  theoretical  value  of  the. coding  error  as  given  by  (1) 

in  terms  of  the  peak-to-peak  signal  to  rms  noise  ratio  for  a  discrete  random 
field  with  an  autocorvariance 


R(x,x,y,y)  =  e"c‘lx-Xl"f3l>’-y| 


(12) 


using  Karhunen-Loeve,  Hadamard,  discrete  Fourier,  discrete  Cosine  [11],  and 
Slant  [12]  transformations  for  c  =  0.0545,  B  =  0.128.  These  curves  are  obtained 
using  a  bank  of  n  DPCM  encoders  with  one-element  predictors.  The  quantizers  in 
the  DPCM  systems  are  of  instaneous-companding  type  and  are  designed  for  the 

*enddresultWrequir1ngmiessEcornputationsC[10].b°th  "  "r  «  «'»•  the  same 
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probability  density  functions  of  the  differential  signals  in  the  DPCM  sys¬ 
tems.  Experimental  results  (see  Figure  4)  indicate  that  the  probability  density 
functions  of  the  differential  signals  are  double-sided  exponential  functions. 

Figure  5  shows  that  increasing  the  block  size  M  improves  the  theoretical 
performance  of  the  proposed  one-dimensional  systems,  however  the  improvement 
becomes  negligible  for  values  of  M  larger  than  8.  This  makes  this  coding  sys¬ 
tem  less  sensitive  to  the  block  size  than  the  standard  transform  coding  systems 
which  use  memoryless  quantizers. 

HI.  Two-Dimensional  Transformation  and  DPCM  Coding 

The  most  severe  limitation  of  the  two-dimensional  transform  coding  tech¬ 
nique  is  the  large  number  of  operations  needed  for  coding  pictorial  data  and 
the  resulting  equipment  complexity.  This  limitation  becomes  less  significant 
when  the  picture  is  divided  into  smaller  block  sizes  then  encoded,  but  this 
also  limits  the  efficiency  of  the  system  since  the  elements  of  various  blocks 
remain  correlated  in  the  transformed  domain.  A  coding  system  utilizing  two- 
dimensional  transformation  and  a  bank  of  DPCM  systems  would  exploit  this  correla¬ 
tion,  thus  improving  the  coding  efficiency  of  the  system. 

In  this  system  the  sampled  image  u(.,.)  of  N  by  N  elements  is  divided  into 
arrays  of  M  by  M  elements.  The  location  of  each  array  in  the  image  is  indexed 
by  subscripts  (M).  A  two-dimensional  unitary  transformation  of  each  block  of 
image  is  obtained  and  is  ordered  to  form  a  one-dimensional  array  U^,  i  =  1,2,..., 
M2  where  Ui (k,£)  refers  the  the  ith  member  of  this  array  which  is  obtained  by 
transforming  the  (k,£)  block  of  the  data.  This  transformation  and  its  inverse 
are  modeled  as 
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(13) 


M  M 

U^(k.A)  =  u(kM+x,£M+y)^(x,y) 

M2 

u(kM+x,£M+y)  =  z  Ui(k,iH.(x,y)  (u) 

where  ♦,(*■»  are  a  set  of  orthonormal  basis  matrices.  The  elements  of  u,(k  i) 

arrays  for  various  values  of  1,  (1-1.2 . M2)  are  correlated,  thus  could  be’ 

coded  by  M  DPCM  ;ystems.  The  analysis  of  this  problem  Is  similar  to  one  dis¬ 
cussed  In  Section  II.  The  coding  error  and  the  bit  assignment  rules  are 

e2  -  R(0.0)  -  £  C.I01  +  fe  ?  .'im1  2  .  . 


W.O)-^  0,(0)  +  £  £  e  1  e* 


Mh  i 

m.  =  —  +  1 
i  n  a 


2  1  n  ? 

£n  e.  -  -  z  in  p 

l  n  ^  xn  ei 


where 


°l(T>  =  y»l  y=i  x=l  J--1  R(x'M  (x-y)*,(S.y)  (17) 

and  e2  is  the  variance  of  the  differential  signal  In  1th  DPCM  system. 

2  m 

ef  '  cj(0)  -  Ai jci (J )  (18) 

Note  that  in  (16),  R(.,.)  1s  the  autocovariance  of  the  Image  u(.,.)  and  It  Is 
assumed  that  u(.,.)  is  a  covariance  stationary  random  field. 

The  last  term  In  (15)  Indicates  the  improvements  due  to  the  use  of  DPCM 
systems.  Replacing  the  DPCM  systems  by  memoryless  quantizers  will  give  a  sim¬ 
ilar  expression  for  e2  where  e2  Is  replaced  by  C,(0).  Thus  the  coding  Improve¬ 
ment  due  to  the  use  of  DPCM  encoders  depends  on  the  relative  values  of  e2  and 


< 

f  ‘ 


c.(0)  which  in  turn  depends  on  R(.,.)  and  the  block  size  M.  In  addition,  note 
that  the  coding  error  of  the  two-dimensional  transform  encoders  are  obtianed  from 
(15)  using  0^(0)  for  «.  and  b  =  1 ,  since  it  is  generally  assumed  that  the  trans¬ 
formed  data  possesses  a  normal  distribution  and  according  to  (8),  b  =  1.  This 
assumption  is  not  accurate  for  small  block  sizes  and  this  partially  accounts  for 
the  difference  between  the  experimental  and  the  theoretical  performance  of  the 
standard  two-dimensional  transform  encoders  [4].  In  the  hybrid  systems  e2  is 
evaluated  for  b  =  2  which  gives  a  larger  theoretical  value  for  the  coding  error, 
however  experimental  results  indicate  that  the  differential  signal  in  the  hybrid 
systems  posses  an  exponential  density  regardless  of  the  block  size  (see  Figure  4). 
This  results  in  a  better  agreement  between  the  theoretical  and  the  experimental 
results  in  the  hybrid  system. 

Figure  6  shows  the  performance  of  the  two-dimensional  system  for  the  ram- 
dom  field  that  was  considered  in  Section  II.  Here  a  block  size  of  4  by  4  is 
used,  the  bit  assignment  is  according  to  (16),  ar.d  the  DPCM  encoders  are  identi¬ 
cal  to  those  employed  for  the  one-dimensional  system.  Figure  7  shows  the  per¬ 
formance  of  the  two-dimensional  systems  for  various  block  sizes  M.  The  results 
are  similar  to  those  observed  for  the  one-dimensional  system. 

IV.  Experimental  Results  and  the  Noisy  Channel 

The  coding  systems  discussed  in  Sections  II  and  III  were  simulated  on  a 
digital  computer  and  were  used  to  encode  the  picture  shown  in  Figure  9a.  This 
picture  is  composed  of  256  by  256  picture  elements,  each  element  quantized  to 
64  levels.  Values  of  a  and  8  for  (12)  were  approximated  from  Figure  9a  and  are 
0.0545  and  0.128  respectively.  The  experimental  results  for  the  hybrid  system 


their  counterparts  encoded  by  the  two-dimensional  transform  coder. 

To  study  the  effect  of  channel  error  in  the  performance  of  the  hybrid  en¬ 
coders  a  binary  symmetric  channel  was  simulated  on  a  digital  computer  to  link 
the  transmitter  with  the  receiver.  This  channel  would  operate  on  each  binary 
digit  .-independently  changing  each  digit  from  0  to  1  or  from  1  to  0  with  probability 
P  and  leaving  the  digit  unchanged  with  probability  1-p.  Then  the  same  receiver 
discussed  in  Section  II  is  used  to  reconstruct  the  encoded  picture  from  the 
string  of  binary  digits  at  the  output  of  the  channel.  Figure  12  shows  the  degrada¬ 
tion  of  the  one-dimensional  hybrid  system  using  Hadamard  transformation  and  the 
DPCM  encoders,  in  terms  of  signa'-to-noise  ratio,  for  bit-error  probabilities  of 

P  =  10  5  p  =  10  5  and  P  =  10"4-  ^  is  observed  that  the  performance  of  the 
hybrid  system  is  affected  very  little  for  a  channel  noise  corresponding  to  a 
bit-error  probability  of  10  4,  however  it  increases  significantly  at  bit-error 
probabilities  of  10  3  and  10  2.  Also,  it  is  noted  that  the  effect  of  channel 
error  is  more  emphasized  at  higher  bit  rates,  particularly  for  bit-error  pro¬ 
babilities  of  10  3  and  10  2.  Figure  13  shows  the  encoded  pictures  using  the 
above  hybrid  system  at  bit  rates  of  one  and  two  bits  per  pixel  for  bit-error 
probabilities  of  10'2,  ICf3,  and  ICf4.  These  figures  show  that  aside  from  the 
exceptionally  high  bit-error  probability  of  10"2,  the  degradation  of  the  encoded 
pictures  because  of  noisy  channel  is  not  significant.  Visible  degradations  are 
in  the  form  of  short  runs  of  16  pixels  (the  block  length  of  one-dimensional 
transformations)  which  is  due  to  the  error  in  the  output  of  the  DPCM  system 
that  encodes  the  lower  components  of  the  transformed  signal  We  note  that  there 

is  no  propagation  of  the  error  in  the  vertical  direction  as  is  commonly  observed 
in  simple  DPCM  systems. 
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The  hybrid  system  we  have  investigated  is  one  that  uses  the  optimum  coeffici- 
entS  V  A2 . An  1n  ttle  predictors  employed  in  both  the  receiver  and  the  trans¬ 

mitter  of  the  system  (see  Figure  1).  These  coefficients  depend  upon  the  statistics 
of  the  particular  picture  that  we  are  encoding.  To  make  the  encoder  Independent 
of  the  signal  statistics,  a  conmon  value  of  unity  can  be  used  for  A,  through  A 
This  makes  the  bank  of  DPCM  systems  similar  to  the  differential  encoders  that 
have  been  used  in  practice  [3],  Figure  14  shows  the  performance  of  the  hybrid 

system  discussed  above  using  a  cornnon  value  of  unity  for  A,,  1=1 . .  at  various 

bit  rates  and  bit-error  probabilities.  This  figure  also  shows  the  performance  of 
the  same  encoder  using  optimum  coefficients  and  a  noiseless  channel.  The  degrada¬ 
tion  due  to  making  these  coefficients  unity  is  aboutonedB  and  It  is  almost  the 
same  at  various  bit  rates.  The  behavior  of  this  system  for  various  bit-error 

probabilities  is  very  similar  to  the  performance  of  the  encoder  using  optimum 
coefficients  shown  in  Figure  12. 

Encoded  pictures  using  this  suboptimum  encoder  for  one  and  two  bits  per 
pixel  at  various  bit-error  probabilities  are  shown  in  Figure  15.  Comparison  of 

Figures  13  and  15  shows  tnat  there  is  no  noticeable  difference  in  the  performance 
of  the  optimum  and  suboptimum  encoder. 

To  compare  the  performance  of  the  hybrid  system  with  the  encoder  using  two- 
dimensional  transform  and  block  quantization  for  a  noisy  channel,  the  binary 
symmertric  channel  discussed  aboveis  simulated  in  the  system  that  uses  two-dimen¬ 
sional  Hadamard  transform  and  a  block  quantizer  for  a  block  size  of  16  by  16. 

Figure  16  shows  the  performance  of  this  system  in  terns  of  slgnal-to-nolse  ratio 
at  various  bit  rates  for  bit-error  probabilities  of  10"4,  10"3,  and  10"2.  Two 
pictures  encoded  by  this  system  using  one  and  two  bits  per  pixel  and  bit-error 
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probability  of  10  3  are  shown  in  Figure  17.  Comparison  of  these  pictures  with 
those  in  Figures  11  and  13  indicates  that  the  subjective  quality  of  the  resul¬ 
tant  pictures  is  affected  more  because  of  channel  noise  when  coded  by  the  hybrid 
system  then  it  is  using  the  two-dimensional  transform  technique,  though  the 
hybrid  system  produces  pictures  with  a  better  overall  quality. 

V.  Conclusions 

The  hybrid  system  using  the  two-dimensional  transformations  and  the  DPCM 
encoders  is  essentially  the  standard  two-dimensional  transform  coding  system 
where  the  encoding  efficiency  is  improved  by  exploiting  the  inter-block  correla¬ 
tion  in  the  data.  The  system  is  attractive  since  a  small  block  size  could  be 

used  reducing  the  computational  complexity  without  affecting  the  performance  of 
the  system. 

The  hybrid  system  that  uses  one-dimensional  transformation  and  DPCM  systems 
attempts  to  uncorrelate  the  data  in  one  spatial  direction  by  the  transform  tech¬ 
nique  and  in  the  other  spatial  direction  by  the  DPCM  encoders.  The  quantization 
and  the  channel  error  are  introduced  in  the  transform  domain,  where  it  is  less 
objectionable  to  a  human  viewer.  The  fact  that  the  system  uses  a  one-dimensional 
transformation,  it  is  not  sensitive  to  the  block  size,  and  allows  parallel  opera¬ 
tion  on  the  data,  makes  the  system  attractive  from  an  implementations!  viewpoint. 
The  significant  result  is  the  improved  performance  of  the  system  as  indicated  by 
both  theoretical  and  experimental  results  that  surpasses  the  performance  of  both 
the  DPCM  and  the  two-dimensional  transform  coding  systems  and  the  fact  that  the 

performance  of  the  encoder  Is  not  impaired  significantly  for  small  or  moderate 
levels  of  channel  noise. 
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cascade  of  one-dimensional  transformations  and 


Bit  rate  versus  the  logarithm  of  the  quantization  error  and  its  linear  approxi¬ 
mation  for  a  variate  with  a  double  sided  exponential  probability  density  functio 
using  an  instaneous -companding  quantizer. 


TWO  D IMENS IONRL 


Theoretical  performance  of  the  propoeed  one-dimensional  encoders  versus 

block  size  M.  The  solid  and  the  dashed  lines  refer  to  the  Karhunen-Loev.  and 
the  Hadamard  transforms,  respectively. 
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SIGNAL  TO  NOISE  RATIO  D8 

;al  results  comparing  the  performance  of  the  proposed  hybrid  systems  using  the  picture 
a).  The  performance  of  a  third-order  simple  DPCM  and  the  encoder  using  two-dimensional 
;ransform  and  a  block  quantizer  is  included  for  comparison. 


(c)  Hadamard 


Figure  10.  The  encoded  pictures  using  the  cascade  of  one¬ 
dimensional  transformations  and  DPCM  systems,  2  bits/ 
pixel,  M  =  16. 
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(a)  1  bit /Pixel 

Two-Dimensional  Hybrid 


(b)  2  bits/Pixel 
Two-Dimensional  Hybrid 
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(c)  1  bit/Pixel 

Two-Dimensional  Hadamard 
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1  'Cure  11.  (a)  and  (b)  are  coded  by  the  hybrid  system  using  two-dimensional 
Hadamard  transform  (4  by  4)  and  DPCM.  (c)  and  (d)  are  coded  by  the  system 
using  two-dimensional  Hadamard  transform  (16  by  1  o)  and  a  block  quantizer. 


B3 


Preceding  page  blank 


UM  HYBRID  SYSTEM 


SIGNAL  TO  NOISE  RATIO  DB 

Performance  (experimental  results)  of  optimum  hybrid  system  (A-j's  have  correct 
one-dimensional  Hadamard  transform  and  DPCM  encoders  for  a  noisy  channel. 


Figure  13.  Coded  by  the  hybrid  system  using  one-dimensional  Hadamard 
transform  and  DPCM  encoders  employing  a  noisy  channel. 
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SIGNRL  TO  NOISE  RATIO  DB 

rig.  14.  Performance  (experimental  results)  of  suboptimum  hybrid  system  (A-j  =  1)  using  one-dimensional 
Hadamard  transform  and  DPCM  encoder  for  a  noisy  channel.  Performance- of  optimum  hybrid 
encoder  is  included  for  comparison. 


Reproduced  from 
best  available  copy. 


ck  quantizer  for  a  noisy  channel. 


(b)  2bits /pixel  ,  P=10 


Figure  17.  Coded  by  the  system  using  two-dimensional  Hadamard  transform 
(16  by  16)  and  a  block  quantizer  with  a  noisy  channel. 
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REAL  TIME  DISCRETE  FOURIER  TRANSFORMS  USING  SURFACE  ACOUSTIC  WAVE  DEVICES 


J.  M.  Alsup,  R.  W.  Means,  and  H.  J.  Whitehouse 
Naval  Undersea  Center,  San  Diego,  California 


Abstract 


A  real  time  Discrete  Fourier  Transform  (DFT)  module  is  presented  using 
surface  acoustic  wave  devices.  The  chirp  -  z  algorithm  is  used  to  compute 
the  DFT.  It  requires  premultiplication  of  the  time  signal  by  a  discrete 
complex  FM  chirp.  Since  the  algorithm  is  complex  it  is  implemented  in 
real  arithmetic  by  four  convolvers  operating  in  parallel.  A  63  tap 
transversal  filter  performs  the  convolution  required  for  a  length  32  DFT. 
Each  transversal  filter  has  uniformly  spaced  taps  whose  weights  are  the 
cosine  or  sine  components  of  a  complex  chirp.  The  computation  time  of 
the  device  is  proportional  to  the  number  of  samples,  N,  in  the  time 
signal.  If  N  is  a  power  of  2  the  computation  time  of  the  Fast  Fourier 
Transform  (FFT)  is  N  log2N.  This  device  is  thus  a  factor  of  log2N  faster 
than  a  FFT  using  a  multiplier  of  the  same  speed.  It  is  also  a  low  power, 
light  weight  device. 
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REAL  TIME  DISCRETE  FOURIER  TRANSFORMS  USING  SURFACE  ACOUSTIC  WAVE  DEVICES 


J.  M.  Alsup,  R.  W.  Means,  and  H.  J.  Whitehouse 


Introduction 

A  transformation  of  considerable  importance  for  signal  processing  is  the 
discrete  Fourier  transform  (DFT).  It  may  be  calculated  using  the  chirp-Z 
transform  (CZT)  algorithm,  [1,2]  a  method  which  is  suitable  for  implemen¬ 
tation  via  acoustic  surface  wave  devices  as  well  as  other  forms  of  the 
transversal  filter.  The  CZT  algorithm  decomposes  the  DFT  into  a  premul¬ 
tiplication  by  a  discrete  chirp,  a  convolution  with  a  discrete  chirp,  and 
a  postmultiplication  by  a  discrete  chirp.  This  is  illustrated  in  Figure  1 
and  by  the  following  equations  for  an  N-point  DFT: 


but  , 


e-j2imm/d 


-2mn  =  -m2  +  (m-n)2  -  n2; 


therefore, 


Since  this  algorithm  requires  the  use  of  complex  numbers,  its  implementa¬ 
tion  with  real  hardware  requires  the  use  of  four  convolvers.  The  real  and 
imaginary  components  of  the  input  g0  are  multiplied  by  the  real  and 
imaginary  components  of  the  premultiplier-chirp  and  combined  in  pairs  to 
drive  tne  inputs  of  four  chirp  convolvers.  The  convolver  outputs  are  com¬ 
bined  in  pairs  and  multiplied  by  the  real  and  imaginary  components  of  the 
postmultiplier-chirp  and  combined  again  to  provide  the  real  and  imaginary 
components  of  the  output  G  . 

The  computation  time  of  the  device  is  proportional  to  the  number  of  samples, 
N,  in  the  time  signal.  If  N  is  a  power  of  2  the  computation  time  of  the 
Fast  Fourier  Transform  (FFT )  is  N  log2N.  This  device  is  thus  a  factor  of 
log2N  faster  than  a  FFT  using  a  multiplier  of  the  same  speed.  It  is  also 
a  low  power,  light  weight  device. 


Transversal  Filter  Convolver 


The  convolution  operation  may  be  implemented  in  highly  parallel  form  as  a 
transversal  filter.  Such  a  filter  is  diagrammed  in  Figure  2  along  with 
the  equivalent  matrix  multiplication.  The  filter  consists  of  a  delay  line 
which  propagates  the  input  signal  gn  along  with  discrete,  equally-spaced 
taps  which  sample  the  signal  along  the  delay  line.  The  taps  are  weighted 
by  the  values  bn  which  correspond  to  the  real  and  imaginary  components  of 


The  authors  are  with  the  Naval  Undersea  Center,  San  Diego,  California 
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the  discrete  chirp:  bn  =  cos  (irn2/N)  or  sin  (irn2/N).  When  N  is  even  the 
tap  weights  repeat  and  circular  convolution  can  be  used.  The  output  of 
the  filter's  delay  line  is  shown  as  connectable  to  its  input  so  as  to  be 
able  to  perform  circular  convolution.  The  complete  DFT  via  the  CZT 
algorithm  is  shown  in  Figure  3.  Parallel  implementation  of  the  complex 
arithmetic  is  shown  for  simplicity. 

Surface  Wave  CZT  Convolver 


For  signal  processing  at  sample  rates  above  5  MHz,  the  acoustic  surface 
wave  technology  is  a  very  suitable  candidate  for  implementation  of  a  CZT 
convolver.  Not  only  can  the  sine  and  cosine  tap  weightings  be  easily 
implemented  (by  adjusting  the  finger  overlap  in  the  interdigital  trans¬ 
ducer  taps),  but  all  four  convolvers  required  for  the  algorithm  can  be 
fabricated  on  a  single  substrate. 

Figures  4  and  5  show  a  CZT  mask  used  at  the  Naval  Undersea  Center  to 
fabricate  a  surface  wave  4-channel  convolver.  In  this  instance,  circular 
convolution  is  simulated  by  using  a  repeating  structure  of  63  taps  to 
implement  a  32-point  transform:  the  last  31  taps  are  a  replication  of  the 
first  31  taps. 

Features  of  the  Design 

The  launch  transducers,  located  at  the  extreme  right  and  left  ends  of  the 
mask,  have  an  aperture  which  is  common  to  two  of  the  tapping  structures. 
Thus,  even  though  there  are  four  receiving  structures  which  constitute  the 
four  transversal  filters,  there  are  but  two  acoustic  paths  each  of  which 
is  occupied  by  a  sine  and  cosine  pair  of  tapping  structures.  Acoustic 
waves  may  be  launched  from  either  end  since  the  tapping  structures  are 
symmetric.  Since  the  launch  aperture  is  about  318  acoustic  wavelengths 
wide  and  the  maximum  propagation  distance  is  about  851  acoustic  wave¬ 
lengths  long,  it  can  be  established  that  the  farthest  tap  is  well  within 
the  near  field  of  the  launch  transducer  and  that  the  acoustic  waves  travel 
across  the  substrate  in  well -colli mated  beams.  The  fractional  'far-field 
length'  involved  is  the  ratio  851 / ( 31 8 ) 2  4  .008. 

The  mask  is  transferred  to  an  ST-cut  quartz  substrate  as  a  pattern  of 
aluminum  fingers  oriented  for  acoustic  propagation  in  the  x-direction. 

The  size  of  the  crystal  is  3.5  inches  by  1.75  inches  by  (1/8)  inch. 

The  launch  transducers  and  each  receive  tap  are  composed  of  three  pai^s  of 
interdigitated  fingers  spaced  at  .0044-inch  intervals.  The  finger  width 
and  the  finger  gap  are  equal  at  .0011  inch.  Since  the  velocity  of  pro¬ 
pagation  is  3158  m/s,  the  resulting  center  frequency  is  28.3  MHz. 

A  close-up  view  of  the  tapping  structure  shows  how  the  finger  groups  are 
weighted  (Figure  5).  The  finger  is  cut  short  of  the  full  possible 
aperture  if  a  weight  smaller  than  one  is  desired.  Finger  tips  are 
retained  to  minimize  possible  acoustic  diffraction  due  to  differential 
mass  loading  across  the  aperture.  A  tap  of  weight  zero  is  implemented  by 
a  space  (no  finger  group). 

The  ratio  of  group  spacing  to  wavelength  was  originally  designed  to  be 
six,  but  artwork  limitations  resulted  in  a  ratio  of  5.91.  Thus,  the 
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inter- tap  delay  became  0.209  psec  and  the  corresponding  chirp  sample 
frequency,  4.78  MHz.  The  resulting  fractional  bandwidth  is  4.78/2(28.3)  = 
.084. 


The  rate  of  calculation  of  the  OFT  is  thus  determined:  32  points  per 
6.69  nsec.  When  calculating  such  a  rate  in  the  manner  described,  it  is 
assumed  that  two  units  are  multiplexed  so  that  one  is  calculating  while 
the  other  is  loading  with  new  data. 

The  ends  of  the  quartz  crystal  were  rounded  so  as  to  permit  propagation 
of  the  acoustic  waves  around  the  ends  of  the  back  side  of  the  crystal 
where  they  are  absorbed  by  a  thin  layer  of  silastic  room- temperature- 
vinyl  (RTV)  distributed  over  the  back  surface. 

The  crystal  was  mounted  in  a  milled  aluminum  box  with  separate  compartments 
in  the  underside  to  accommodate  electrical  matching  circuits  for  each 
input  and  output  signal  with  large  RF  isolation.  The  coaxial  connectors 
are  SMA  female.  The  top-side  was  fitted  with  a  transparent  plastic  top 
inside  the  standard  metal  top  so  as  to  afford  more  protection  under 
routine  inspection.  Leads  were  bonded  from  feedthrough  terminals  to  the 
surface  electrode  busses  via  conductive  silver  paint. 

Preliminary  Experimental  Results 

Two  devices  were  constructed.  One  was  used  to  generate  the  32-point  sine 
and  cosine  chirp  samples  for  premultiplication,  and  the  other  to  convolve 
the  multiplied  input  with  the  four  63-point  chirp  transversal  filters. 

Values  used  to  construct  the  weights  in  the  tapping  structure  are  given 
in  the  first  three  columns  of  Table  1.  The  impulse  response  of  a  32-point 
tapping  structure  is  shown  in  Figure  6.  A  bulk-wave  cancellation  scheme 
was  used  to  obtain  this  desired  output,  since  some  spurious  signals  were 
otherwise  obtained  with  the  electrode  pattern  on  ST-cut  quartz. 

Values  expected  at  the  output  terminals  of  the  appropriately-sunned  con¬ 
volver  electrodes  are  given  in  the  last  three  columns  of  Table  1  for 
the  case  when  the  input  g  is  a  constant  =  1.  This  case  simplifies  the 
evaluation  since  two  of  the  convolver  operations  and  the  entire  post- 
multiplier  operation  become  unnecessary.  Values  actually  observed  (after 
the  loading  interval)  are  shown  in  Figure  7,  and  the  final  summed  output 
in  Figure  8. 

It  is  expected  that  the  next  generation  device  will  have  a  better  time- 
domain  response  such  that  the  desired  cancellation  of  the  component 
responses  will  be  more  nearly  perfect  in  the  region  following  the  dc  peak. 

Concl usions 


The  acoustic  surface  wave  implementation  of  CZT  convolvers  and  discrete 
chirp  generators  appears  to  hold  much  promise.  Though  the  outputs 
observed  to  date  are  not  yet  satisfactory,  they  are  expected  to  improve 
w"'th  refinement  and  practice.  A  similar  surface  wave  generator/convolver 
utilizing  nearly  identical  design  techniques  has  been  shown  to  perform 
with  a  very  high  level  of  accuracy  [31. 


ofHsamDlP^T^^jr  r-me  °-  th?  algo!:ithm  is  Proportional  to  the  number 
o  samples,  N,  in  the  time  signal,  a  time-domain  implementation  such  as 

YnVJ.eXp?cted  *?  speed  up  calculation  of  the  OFT  by  a  facto?  of 
log2N  relative  to  an  FFT  with  identical  capacity  and  speed. 

The  surface  wave  implementation  is  relatively  small  and  light-weight,  and 

at  Mnh»  Jequlrem®nt  forA  electrical  power.  The  potentials  for  operating 
at  higher  frequencies  and  over  longer  acoustic  paths  (as  compared  to  the 
design  parameters  reported  here)  point  toward  larger  capacity,  faster 
acoustic  surface  wave  CZT  devices  in  the  future. 
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Fig.  1.  Chirp  Z  Transform  Implementation  of  the  DFT 
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Table  1.  32-Point  CZT:  Filter  Tap  Weights  &  Convolution  Output  Values  for  DC  Input 


FIGURE  6.  SURFACE  WAVE  CZT  FILTER  32-POINT 
IMPULSE  RESPONSE:  (a)  sin  7rnz/32;  (b)  cos  7rn2/32. 
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MOURE  7,  SURFACE  WAVE  CZT  OUTPUT  COM 
PONENTS  FOR  DC  INPUT:  (a)  sin(32)  *  sin  (63); 
(b)  cos(32)*cos  (63). 


FIGURE  8.  SURFACE  WAVE  CZT  OUTPUT  FOR 
DC  INPUT. 
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ABSTRACT 

A  real  time  Discrete  Fourier  Transform  (OFT)  device  is  presented  using  bucket 
brigade  devices.  The  chirp  -  2  algorithm  is  used  to  compute  the  DFT.  It  req¬ 
uired  premultiplication  of  the  time  signal  by  a  complex  FM  chirp,  convolution 
with  a  complex  FM  chirp  and  postmultiplication  by  a  complex  FM  chirp.  Since 
the  algorithm  is  complex  it  is  implemented  in  real  arithmetic  by  four  convolvers 
operating  in  parallel.  A  two  hundred  stage  bucket  brigade  transversal  filter 
performs  the  convolution  required.  The  signal  is  tapped  with  a  source  follower 
Whose  load  determines  the  weighting  coefficient.  The  computation  time  of  the 
device  is  proportional  0  the  number  of  samples,  N,  In  the  time  signal.  The 
computation  time  of  the  Fast  Fourier  Transform  (FFT)  is  N  log2N  if  N  is  a  power 

Of  2.  The  device  is  thus  a  factor  of  log2N  faster  than  the  FFT.  It  is  also  a 
low  power,  low  weight  device. 


I.  INTRODUCTION 


Because  of  their  analog  nature,  charge  coupled  devices  (CCD's)  offer  the  pos¬ 
sibility  of  performing  many  of  the  signal  processing  functions  which  until  now 
have  fallen  into  the  realm  of  digital  filtering.*  The  discrete  Fourier  trans¬ 
form  (DFT)  of  an  electrical  signal  is  one  such  linear  filtering  operation  which 

readily  lends  itself  to  implementation  with  CCD's  and  is  discussed  in  this 
paper. 


e  technique  of  using  CCD's  to  perform  the  discrete  Fourier  transform  (DFT)  of 

an  electrical  signal  was  first  proposed  by  Whitehouse,  Speiser,  and  Means  2 

Ihe  technique  they  have  proposed  is  called  the  chirp-Z  transform  (CZT)  because 

tne  mathematical  operations  involved  multiplication  and  convolution  with  chirp 
signals.3 


The  CZT  was  first  conceived  as  an  algorithm  for  the  DFT  on  a  digital  computer 
It  never  gained  wide  acceptance,  however,  because  the  number  of  digital  oper¬ 
ations  required  for  its  implementation  is  the  same  as  for  the  Cooley-Tukey  fast 
Fourier  transform  (FFT)  algorithm  It  can  be  viewed  as  three  distinct  ooerations 
all  of  which  can  be  perfumed  with  analog  circuitry.  (1)  prmultiplication  by 
a  chirp  waveform,  (2)  convolution  with  a  chirp  waveform,  and  (3)  post-multipli- 

cation  by  a  chirp  waveform.  These  three  operations  are  illustrated  in  the  block 
diagram  of  Figure  1 . 

The  DFT  of  a  band  limited  electrical  signal  g(t)  is  found  by  sampling  the  signal 

at  uniformly  spaced  instants  of  time  to  obtain  the  N  points  g  .  The  DFT  is 
then  given  by  the  definition  n 


N-l 

6k  =  * 
n=o 


gn  exp[- 


c  <  k  <  N-l 


(1) 


The  CZT  algorithm  for  evaluation  equation  (1)  is  derived  by  making  the  substitu- 
tion 


2kn  =  k?  +  n2  -  (n-k)2 
into  equation  (1)  to  get 


(2) 


95 


(3) 


Gk  s  e-lTTk2/N 


N-l 


I 

nBo 


e t-r(n-k) 2/N  e-lTrn2/N 

9 « 


quat,0"  (3>  Su«ests  Owl*  the  hardware  fomentation  Illustrated  tn 

Ze  sal  Hltal9Tthm  °ff6rS  U"1qUe  adMnt89eS  Whe"  1mp“  using  CCD 
Z  5  eCaUSe  *"  "*  Pr°CeSS1"9  “»  *  performed  ,n  an  an  log 


:::l:  "*  T  ?er-9lVeS  eXpertaM1  verification  of  the  CZT 

c  n  ept.  The  der,onstrat,on  Is  only  approximate  but  It  Illustrates  the  simplicity 

IV  concludTwUh  the4""  ''V65"1665  ^  ““ct  C2T  ^f'^ntatlons.  Section 
concludes  with  the  discussion  of  applications. 


XL — preliminary  expfrimental  results 


To  demonstrate  the  concept  of  the  CZT,  the  simplified  system  shown  at  the  too 
igure  was  Implemented  for  approximating  the  power  density  spectrum  The 

the  functions!"  ^  de“"Strati0"  h"»  -ponses  which  approximate 


h(t)  =  cos  ut2 
for  the  COS  filter  and 


Td.  Td 
—  <  t  <  — 


(4) 


g(t)  =  sinpt2 


<  t  < 


2 


(5) 


dpi  flUer‘  They  are  bucket  bri9ade  device  ( BBD )  filters  which  were 

e  e  oped  for  a  spread-spectrum  recejver .v,s,e,z  The  impulse 

re'at,0n  °f  ‘*>ese  filters  are  given  in  Figure  3. 


In  the  deration,  the  chirp  signal  shown  in  Figure  2  chirps  from  97  5  kHz  to 

102.5  kHz  in  5  msec,  and  gives  the  approximate  power  density  spectrum  !  ‘ 

frpQupnrv  ka tu  , .  *  spsctrum  3CPOSS  this 

way.  If  the  in  e.°Pera  10n  of  the  system  ca"  oe  understood  in  the  following 

..  '  t'  PU  S19na  ls  a  single  frequency  sinusoid  of  frequency  f  then 
the  output  of  th»  mixpr  hac  a  w-ir-p  equency  r  ,  then 

to  f  102  c  k«  '  difference  component  which  chirps  from  f  -97  5  kHz 

s  102.5  kHz  and  a  sum  component  which  is  filtered  out  in  a  low-pas's  filter 

after  mixing  (not  shown  In  Fio„re  21  n..  .  .  pass  filter 


after  mixing  (not  shown  In  Figure  2)  T  e  f  ,  ^  f1Uer 

ngure  2).  The  filters  are  clocked  at  20  kHz  and  have 
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impulse  responses  which  chirp  from  -5  kHz  to  +5  kHz  in  10  msec, 


The  convolution  operation  is  indicated  schematically  at  the  botton  of  Figure  ■> 
The  frequency  of  the  filter  impulse  response  is  indicated  by  the  line  extending 
rom  -5  kHz  to  +5  kHz,  The  input  signals  are  reversed  in  time  as  required  by  ' 

T nP  rnrn/nl  1 14-  -J  j. .• _  .  1  J 


,  ,  - - ■  "  ‘'""c  uo  i  c^u  i  r  cu  Dy 

the  convolution  operation,  and  are  shown  for  four  values  of  f  equally  spaced 


w 1  1  -  cuuaiiy  bUdCea 

across  the  chirp  band:  97,5  kHz,  99.17  kHz,  100.63  kHz  and  102.5  kHz.  A  cor¬ 
relation  peak  in  the  filter  output  occurs  when  the  frequency  of  the  time  reversed 
input  signal  "lines  up"  with  the  frequency  of  the  impulse  response.  At  the  1 
instant  shown,  the  97.5  kHz  signal  lines  up  with  the  filter  and  gives  a  cor- 

la6  Ir'in  "T'  "'****  ^  COrre,ati°"  P«k«  which  occur 

later  in  time  by  the  amount  of  At  given  by 


At  =  £(f  -97.5kHz) 


(6) 


Whe<  is  the  slope  of  the  frequency  vs  time  curve  as  can  he  seen  from  egs.  (4) 


Two  filters  are  used  to  perform  in-phase  and  quadrature  processing,  and  the 
output  signal  which  results  from  squaring  and  summing  the  filter  outputs  closely 

r™kHz  ha”  denS''ty  SP6Ct™  °f  1nPUt  Si9"a’  aCr°SS  «•  97-9  ^z 


Z7?-"  R9:r:  2  d°eS  "0t  9iVe  the  exact  P0Wer  de"5t*  =PP«rum.  There  are 
purious  sidelobes  due  to  an  approximation  made  that  an  up  chirp  is  ortho 

^input's  Ii^tT-  I"9  °UtPUt  haS  3  dePendenCe  UP°"  the  Phasi-  of 

t  me  ha  d  a  P  ""  ""  be  redUCed  b*  ™ters  »f  larger 

-bandwidth  product  or  by  implementing  the  exact  in-phase  and  quadrature 
processing  to  be  discussed  in  Section  III.  quadrature 


The  actual  operation  of  the  system  is  Illustrated  in  the  series  of  photographs 
shown  ,n  Figure  4.  The  input  is  a  single  frequency  sinusoid  whose  frequency 

indkatad  in  - «. outPUt  shows  a  a[irr 

correlation  peak  whose  position  scales  uniformly  with  input  frequency  The 
time  axis  calibration  is  1.67  msec/cm,  and  the  correlation  peak  shifts  hv  1 
each  time  the  input  frequency  increases  by  1.67  kHz  („/,  =  i0e  sec-Z)  ?he  Cm 
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ft. 


correlation  peaks  should  all  be  of  equal  amplitude.  However,  because  of  the 
above  mentioned  dependence  on  input  phase,  there  is  a  small  random  variation 
from  one  photograph  to  the  next. 


III.  FOURIER  TRANSFORM  ANALYSIS 

The  exact  DFT  can  be  implemented  by  the  hardware  shown  in  Figure  5.  The  com¬ 
plex  chirp  Z  algorithm  as  defined  In  equation  (3)  is  implemented  by  real  hard¬ 
ware  operating  in  parallel.  The  transversal  filters  discussed  in  Section  II 
can  implement  a  cosine  transform8  which  is  a  subset  of  the  DFT.  In  order  to 

illustrate  this  it  is  first  necessary  to  discuss  the  symmetries  involved  in 
Fourier  transforms. 


If  the  input  time  signal  g(t)  is  real,  as  it  always  is  in  signal  processing 
applications,  then  the  output  Fourier  transforms  have  the  property  that 


This  is  an  awkward  symmetry  to  make  use  of  in  reducing  data  handling  requirements. 

However,  if  the  input  data  is  real  and  symmetric,  then  the  output  Fourier  trans- 
form  is  also  real  and  symmetric. 

This  symmetric  transform  is  sometimes  referred  to  in  the  literature  as  the  cosine 
transform.  It  is  not  uniquely  defined  and  different  authors  have  defined  slightly 
different  cosine  transforms.  There  are  two  implicit  symmetrizations  possible  on 
the  input  data:  (1)  The  data  set  of  length  N  is  reflected  about  a  data  point. 

This  will  result  in  a  data  set  of  length  2N-1;  (2)  The  input  data  set  of  length 

N  is  reflected  so  that  the  resultant  data  set  is  of  length  2N.  An  example 
shoulo  clarify  these  two  casses. 

1)  Original  data  set:  [1,2, 3, 3, 4] 

Symmetrized  data  set:  [4, 3, 3, 2,1 ,2, 3, 3, 4] 

2)  Original  data  set:  [1,2, 3, 3, 4] 

Symmetrized  data  set:  [4, 3, 3, 2, 1 ,1 ,2, 3, 3, 4] 
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The  Naval  Undersea  Center  is  interested  in  implementing  the  cosine  transform 
because  it  is  Involved  in  a  prog i am  of  television  band  width  reduction.  The 
ordinary  discrete  Fourier  transform  has  several  faults.  Since  it  is  really  a 
Fourier  transform  of  a  periodically  repeated  signal,  it  may  have  artificial 
discontinuities  at  the  boundaries.  If  the  signal  is  made  symmetric  then  the 
discontinuities  disappear  and  the  Fourier  transiorm  will  thus  converge  more 
rapidly.  It  will  be  shown  later  that  making  the  signal  symmetric  does  not 

take  any  hardware.  In  fact  it  reduces  the  amount  of  hardware  required  to 
perform  a  Fourier  Transform. 

In  order  to  demonstrate  this  it  is  advantageous  to  go  through  several  stages 

of  analysis.  Consider  a  continuous  Fourier  transform  of  a  finite  time  signal 
g(t) 


G(w)  =  J*  g(t)e"lwt  dt 


(8) 


Define  the  symmetrized  time  signal  as 


q(t)  =  9(t)  0  <  t  <  T 

9(tJ  -  g(-t)  -T~<  T<  0 


(9) 


Then  the  Fourier  transform  of  g(t)  is 
T  „ 

G(w)  =  j  g (t )  cosuit  dt 
-T 

Since  g(t)  is  real  and  symnetric  then  G(ui)  is  real  and 
T  . 

G(w)  =J  g (t )  coswt  dt 


(10) 


symmetric  and  thus 


-T 


=  2/  g (t)  cosuit 


dt 


OD 


For  future  comparison  it  is  convenient  to  write  thi 


s  as 


G(w)  =  2  Re  f  g(t)e"1'jt  dt 


(12) 


Notice  that  the  Fourier  transform  of  the  symmetric  time  signal  g‘(t)  has  been 
derived  by  only  working  with  the  non-symmetric  time  function  g(t)  and  that 
only  an  integral  of  length  T  was  required. 
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Let  us  consider  now  a  sampled  signal  gn  of  N  points  and  Us  distrain  Fourier 
i-  ansform.  Define  the  svmmotri/'  ?,,***■< 


1  'ans'°nn.  Define  the  symmetric  function 

9nSgn’1  n~  “N+1  <  n  <  N  (13) 

9-n  n^o 

Thns  function  can  be  looked  on  as  symmetric  about  n  =  1/2.  Its  "Fourier" 
form  is  thus 


trans- 


p  N  *  ~2nik(n-l/2) 

°k  =nLNtl  9n  e  35  -N+1  i  k  <  N  (14) 

A 

Since  g  is  symmetric,  then  Im  G.  =  0  and  G  =  G  and  it  ic 

compute^alf  the  set.  k  k  §k+1  $  0nly  "eCeSSary  t0 


;  N  -v  -2ffik(n-1/21 

Gk  =  E  9n  e  2lr“ 

K  n=-N+ln 


o  <  k  <  N-l 


By  rearranging  the  summing  order  and  using  the  fact  that  g  =  g  for  n  >  1 
this  can  be  shown  to  be  n  n-l  -  ’ 


=  2  £  g  cos 

A-A  '  ' 


This  can  be  written  in  the  form 


-2Trk(n+l/2 


o  <  k  <  N-l 


G.  =  2  Re  z  g 


^^Triktntl 
n  2N 


o  <  k  <  N-l 


Notice  that  only  the  original  sampled  time  function  {g  ,  is  used  and  that  the 

urn  is  on  y  over  N  points.  Equation  17  can  be  implemented  by  a  chirp  Z  trans¬ 
form  by  substituting  the  identity. 

2k(n+l/2)  =  k2  +  (n+1/2)2  - (k- (n+1/2) ) 2  (,8) 


This  gives 


-ink2  N-l  J  ,  , 

Gk  =  2  Re  eUT  z  e+l7r^k^n+1/2))2/2N 
n=o 

-i7r(n+l/2)2/2N 


This  cosine  transform  can  be  implemented  by  the  bucket  brigade  devices  by  letting 
N  =  199  except  for  a  factor  of  2  in  the  total  length  of  the  filter.  This  implies 
that  only  the  first  100  cosine  transform  coefficients  will  te  computed  and  that 
the  data  is  defined  to  be  equal  to  zero  after  n  *  100. 

IV.  APPLICATIONS 

The  CZT  algorithm  is  conceptually  a  real  time  algorithm  for  a  sampled  data 
sytem.  The  set  of  points  {gn >  can  be  considered  as  a  time  sampled  function  of 
a  continuous  function  g{t)  of  length  T.  The  transversal  filter  of  the  chirp  Z 
transform  will  have  a  basic  clock  frequency  equal  to  the  sampling  frequency.  It 
has  2N  (or  2N-1)  stages  so  it  is  approximately  twice  as  long  as  the  input  signal. 
A  time  T  must  pass  before  the  set  gp  fills  the  first  half  of  the  transversal 
filter.  Thereafter,  one  DFT  coefficient  is  computed  for  every  time  step.  This 
is  what  is  meant  by  a  real  time  device,  i.e.  for  a  data  sequence  of  length  N 
it  takes  N  clock  periods  to  compute  the  DFT. 


By  this  definition  the  chirp  Z  algorithm  can  be  compared  with  other  algorithms. 
The  basic  algorithm  as  defined  in  equation  (1)  takes  N2  steps  to  compute.  If  N 
is  a  power  of  2,  then  the  Cooley-Tukey  algorithm  called  the  Fast  Fourier  Trans¬ 
form  (FFT)  required  N  log2  N  steps  to  compute.  The  transversal  filter  implemen¬ 
tation  of  the  chirp  Z  algorithm  achieves  its  high  speed  not  by  reducing  the 
number  of  multiples  like  the  FFT  but  by  putting  the  equivalent  of  N  parallel 
multipliers  all  on  the  same  chip.  This  reduces  the  number  of  computational 
steps  from  N2  to  N. 

A  real  time  Fourier  transform  device  will  find  many  applications  in  radar,  sonar 
and  pattern  recognition.  One  of  the  most  promising  is  the  bandwidth  reduction 
of  television.9 
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Fig.  1.  Chirp  Z-transform  Implementation  of  the  DFT 


07.5  -*102.5  kHz 


Fig.  2.  top:  Block  D  agrain  ot  a  Configuration  for  Approx¬ 
imating  the  Power  Density  Sprectrum.  The  filters  marked 
COS  and  SIN  are  clocked  at  20  kHz  so  that  their  impulse 
responses  chirp  from  -5  kHz  to  +5  kHz.  A  low  pass  filter 
v/hich  is  not  shown  follows  the  mixing  operation  and 
eliminates  the  sum  frequency  components. 

BOTTOM:  Schematic  of  the  Convolution  Operation  Performed  • 
Within  the  Filters.  (See  text). 
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Algorithm  with  Parallel  Implementation  of  Complex  Arithmetic 
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APPLICATION  OF  IMAGE  EXTRAPOLATION 
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Abstract 


A  generalized  technique  is  developed  for  estimating  the  missing 
spectral  components  of  the  mathematic?!  transformation  of  «  data  vec¬ 
tor.  The  method  has  application  in  transform  data  coding  systems  and 
time  series  analysis  problems. 

The  spectrum  extrapolation  matrix  operator  is  derived  upon 
principles  of  minimum  mean  square  error  estimation.  Performance  of 
the  extrapolation  method  is  analyzed,  and  experimental  results  of  the 
application  of  the  method  to  one  and  two  dimensional  signals  is  presented. 


Ill 


1.  Introduction 

In  transform  data  coding  systems  one  of  the  coding  techniques  is  to 
transmit  only  those  data  transform  coefficient!  of  low  frequency  or  those 
components  of  low  magnitude  [l  1  At  the  decoder  it  has  been  common  prac¬ 
tice  to  insert  zero^  for  the  missing  components  prior  to  the  inverse 
transformation  operation  that  reconstructs  the  data.  However,  it  is  known 
that  for  all  transformations,  other  than  the  Karhunen-Loeve  transform, 
the  coefficients  are  partially  correlated.  Hence,  it  should  be  possible  to 
utilize  knowledge  of  til  is  correlation  to  estimate  the  missing  coefficients 
rather  than  arbitrarily  setting  their  values  to  zero.  Another  relate'1  pro¬ 
blem  occurs  in  time  series  analysis  in  which  it  is  often  required  to  esti¬ 
mate  the  power  spectral  density  of  a  process  given  only  a  few  points  of 
the  autocorrelation  function  [2].  A  general  solution  to  this  class  of  pro¬ 
blems  has  been  formulated  within  the  framework  of  linear  minimum  mean 
square  estimation  theory. 

2.  Generalized  Spectra 

Consider  an  N  element  vector  i_  which  is  a  sample  of  a  random  pro¬ 
cess  with  zero  mean  and  known  covariance  matrix  K,.  This  vector  under- 

— f 

goes  a  linear  transformation  defined  by  the  N  x  N  element  operator  matrix 
A.  Typical  transformations  of  interest  include  the  Fourier,  Cosine,  Hada- 
mard,  Haar,  Slant,  and  Karhunen-Loeve  transforms  [3-8].  The  N  element 
vector  £  denotes  the  one  dimensional  transformation  of  f_  as  given  by 

£  =  A  f  (1) 

where  the  elements  of_£  represent  the  spectral  components  of  an  expansion 
of  £  in  terms  of  a  set  of  basis  vectors  which  are  the  rows  of  A.  Since  £  is 
zero  mean,  also  possess  a  zero  mean.  Also,  for  a  unitary  transformation 
it  is  easily  seen  that  the  covariance  matrix  of  the  transformed  vector  is 

Kj.  =  A  Kj  A"1  (2) 
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3.  Spectrum  Extrapolation-Noise  Free  Condition 

Figure  la  contains  a  diagram  illustrating  the  generation  of  the  gen¬ 
eralized  spectral  components  of  a  signal,  the  selection  or  truncation  of 
these  components,  and  the  estimation  of  the  discarded  components.  In  the 
diagram  the  complete  vector  of  spectral  components  £  is  operated  upon  by 
an  M  x  N  ( M<  N)  selection  matrix_S  that  extracts  certain  elements  of  £  and 
records  them  in  an  M  element  vector  .  The  selection  matrix  contains 
a  unit  element  at  S(i,  j)  to  transfer  the  jth  element  of  £  to  the  i  th  element  of 
This  selection  matrix  represents  the  sample  selector  of  a  transform 
data  coding  system,  or  it  can  be  considered  to  represent  the  incomplete 
knowledge  of  a  correlation  function  in  the  problem  of  power  spectral  density 
estimation. 

In  the  diagram  of  figure  la,  the  vector  £  is  multiplied  by  an  N  x  M 
estimation  matrix  W  which  provides  an  estimate,  £,  of  the  complete  set  of 
spectral  components  of  £  .  An  inverse  unitary  transformation  A  *  then  re¬ 
constructs  an  estimate,  f  ,  of  the  signal  vector,  {_.  In  essence,  the  esti¬ 
mator  linearly  extrapolates  the  values  of  the  observation  of  known  spectral 
components  to  determine  the  unknown  spectral  values.  The  extrapolation 
is  meaningful  because  of  the  information  provided  by  the  expected  correla¬ 
tion  of  the  elements  of  the  observation. 

The  estimator  matrix  Wis  chosen  to  minunize  the  mean  square  error, 
e,  between  the  signal  vector  and  its  estimate  as  defined  by* 

e  =  trCE((f-f)(f-f)  }]  (3) 

The  optimum  choice  of  Wto  minimize  the  mean  square  error  can  be  found 

A 

by  forcing  the  dynamic  error  (f_-f)  to  be  orthogonal  to  the  observation  jf  [9]. 

notation:  E{  .  }  =  expected  value;  tr=trace  of  a  matrix;  *  =  complex  conjugate; 

T  = transpose, 
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Thus,  by  setting 


E  {  }  =  0  (4) 

where  0_  is  an  NxM  matrix  of  zeros,  one  obtains 

E  C  if-A”1  w  S  A  f)  (S  A_f)*T  3  =  0  (5) 

This  leads  directly  to  the  optimum  solution 

W  =  A  Kf  A"1  ST  fe  A  K  A"1  ST  ]_1  (6) 

Pre multiplying  W  by  the  selection  matrix  S  yields 

S  W  =  iM  (7, 

where  1^  denotes  an  MxM  identity  matrix.  Since 

A 

1  =  w£r  (8) 

it  is  clear  from  equation  (6)  that 

A 

Si  *  S  wjj,  =iT  (9) 

Thus,  the  estimator  simply  copies  the  elements  of  j_  into  the  appropriate 
*  ^  T 

elements  of  J.  The  remaining  elements  of  £  are  obtained  by  a  linear  com¬ 
bination  of  the  elements  of  Examination  of  the  extrapolation  matrix 
computed  from  equation  (6)  indicates  that  the  first  M  rows  of  Wform  an 
identity  matrix.  The  remaining  N-M  rows  of  W  contain  elements  whose 
weightings  are  proportional  to  the  correlation  between  the  coefficient  to 
be  estimated  and  the  spectral  observations.  For  a  Karhunen-Loeve  trans¬ 
formation,  all  transform  coefficients  are  uncorrelated,  and  hence  spectrum 
extrapolation  is  not  applicable.  In  this  instance  the  bottom  N-M  rows  of  W 
are  all  composed  of  zero  elements. 
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The  spectral  estimation  problem  diagrammed  in  figure  la  could 
have  been  formulated  initially  to  estimate  the  signal  f  directly  from  the 
spectral  observation  Let  this  estimator  be  denoted  as  U  where 

t  -  Ur  (10) 

However,  since  A  is  a  unitary  matrix,  the  measure  of  the  error  is  pre¬ 
served  by  the  transformation,  and  it  is  easy  to  show  that  ' 

2  -  A'*  S  (ii) 

Thus,  the  original  solution  is  optimum. 


4.  Spectrum  Extrapolation  -  Additive  Noise 

In  physical  situations  the  observation  usually  is  contaminated  by 
noise.  If  the  noise  is  additive,  as  in  figure  lb,  the  optimum  estimator  that 
minimizes  the  mean  square  error,  c  ,  will,  in  effect,  perform  an  additional 
ask  of  noise  filtering.  Applying  the  orthogonality  principle,  one  obtains 


E  t  (l*i)(v)*T3  =  £ 


(12a) 


E  {  [f-A-1  W  (n+S  A  f)]  [„*T  +  f*T  A'1/  ]  )  =  0 


(12b) 


For  the  usual  case  of  zero  mean  noise  in  which  the  signal  and  noise  are 
correlated,  the  optimum  estimator  is  found  to  be 


W  =  A  Kf  A"1  ST  [£  A  K  A"1  ST  +  K  j'1 


where  K  is  the  noise  covariance  matrix. 
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5.  Computational  Considerations 

In  the  noise  free  case  of  spectrum  extrapolation,  it  is  possible  to 
reformulate  the  expression  for  the  extrapolation  matrix  of  equation  (6)  to 
permit  computation  usin^  a  pseudoinverse  ClO]  algorithm.  Let 

B  s  SAH{  (14) 

where  H,  is  obtained  by  a  factorization  of  the  covariance  matrix  K  .  That 
1  — f 

is 

Kf  ■  HfHfT  (15) 

A  general  factorization  technique  is  given  by 

^  *  easm  (16) 

where  is  a  matrix  containing  the  eigenvectors  of  K^as  column  vectors, 
in  a  diagonal  matrix  whose  elements  are  the  square  roots  of  the  eigen¬ 
values  of _K^ ,  and  M  is  an  arbitrary  unitary  matrix,  which  may  be  chosen 
to  force  to  be  of  lower  triangular  form.  With  this  factorization  the  ex¬ 
trapolation  matrix  expression  becomes 

W  =  A  Hf  BT  (BBV1  (17) 

where  it  is  observed  that  the  matrix  P  =  13  (I3BA)~  can  be  computed  as  the 
pseudoinverse  of  13.  Algorithms  exist  for  the  efficient  computation  of  a 
matrix  pseudoin verse  (jll]. 

6.  Two  Dimensional  Formulation 

The  spectrum  extrapolation  technique  can  be  extended  to  two  dimen¬ 
sional  transforms  quite  easily.  Consider  a  PxP  data  matrix  F  which  is  column 
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scanned  to  obtain  the  N=P  data  vector  f_.  Substitution  may  the/:  be  made 
directly  into  equation  (1),  to  apply  the  subsequent  one  dimensional  process¬ 
ing  algorithm.  A  special  situation  of  interest  occurs  when  the  data  covar¬ 
iance  matrix,  transform  matrix,  and  selection  matrix  are  all  of  separable 
form  such  that  the  following  direct  product  relations  hold 


K.  =  Kr  xK. 
—  f  —  fx  ~fy 


(18a) 


=  A  x  A 
-x  — y 


(18b) 


=  S  x  S 
-x  -y 


(18c) 


where  K,  and  K,  are  (PxP)  row  and  column  covariance  matrices,  A 

—fX  — xV  '  —} 


and  A  are  row  ana  column  transform  matrices,  and  S  and  S  are  row 
-y  -x  -y 


and  column  selection  matrices.  For  the  noise-free  case  it  is  eaaih 


shov/u  tha» 


W  =  W  XW 
—  -x  -y 


(19) 


where 


W  =  A  K  A"1  ST  CS  A  K  A”1  ST] 

-  X  “  X  “  fx  “  X  —x  “X  —  X  —  fx  ~  X  — X 


W  =  A  K  A_1STfe  A  K  A'1  s’1]"1 

-y  —  y  —  fy  —  y  —  y  'Ey  —  y  —  fy  —  y  — y 


In  this  case  spectrum  extrapolation  can  be  sequentially  performed  on  rows 
and  columns  of  the  truncated  two  dimensional  array  of  transform  coeffi¬ 
cients,  thus  permitting  a  significant  reduction  in  computational  complexity. 
Unfortunately,  the  two  dimensional  extrapolation  matrix  cf  equation  (13) 
for  a  noisy  observation  is  not  separable  when  the  separability  conditions 
of  equation  (18)  hold,  even  if  the  noise  is  separable  and  uncorrel&ted. 
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7.  Performance  Evaluation 

The  usual  alternative  to  spectrum  extrapolation  is  simply  to  set  the 

values  of  the  discarded  spectral  components  to  zero  and  perform  an  inverse 

transformation.  For  svch  a  strategy  the  reconstruction  matrix  is  simply 
T 

W  =  S  and  the  estimated  vector  becomes 

1  =  A"1  ST  S  A_f  (20) 

Substitution  into  equation  (3)  yields  the  sub- optimal  mean  square  error 

■*„„  »  -A‘1ST  SA)]  (21a) 

or 

*BO  ■  f  STS)]  (21b) 

If  the  optimal  reconstruction  matrix  of  equation  (6)  or  equation  (13)  is  em¬ 
ployed  the  resulting  mean  square  error  is  found  to  be 

(22) 

Therefore,  a  reduction  in  mean  square  error  by  the  amount 

*  *B0  -  -  ST)  S]  (23) 

can  be  realized  by  spectrum  extrapolation. 

Figure  2  contains  plots  of  the  mean  square  error  of  tx-ansform  com¬ 
ponents  with  and  without  spectrum  extrapolation  for  several  transforms. 

In  this  example  the  data  vector  is  of  length  32  and  the  last  24  elements 
have  been  truncated.  The  data  is  modelled  as  a  first  order  Markov  pro¬ 
cess  with  an  adjacent  element  correlation  of  0.  9”.  The  figures  indicate 
significant  reductions  in  mean  square  error  are  obtained  for  the  Hadamard, 
Haar  and  Slant  transforms.  No  error  reduction  is  possible  with  the 
Karhunen-Loeve  transform  since  the  transform  coefficients  are  uncorre¬ 
lated,  and  there  is  no  information  available  for  extrapolation.  There 
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is  almost  no  improvement  obtained  with  the  Cosine  transform  because  of 

its  close  similarity  to  the  Karhunen-Loev.  transform  for  Markov  process 
data. 


8.  Experimental  Results 

Computer  simulation  experiments  have  been  performed  to  evaluate 
the  spectrum  extrapolation  process.  Figure  3  contains  an  example  of  spec 
tram  extrapolation  in  the  Hadamard  transform  domain  for  two  one-dimen¬ 
sional  signals:  a  sin  tfe  waveform  and  a  pulse  superimposed  on  a  quadratic 
background.  Figure  4  illustrates  two  dimensional  spectrum  extrapolation  h 
the  Hadamard  domain.  In  this  example  a  256x256  pixel  image  was  trans¬ 
formed  in  16x16  pixel  bocks.  The  transform  coefficients  lying  outside  of 
8x8  and  5x5  low  sequency  zones,  for  the  two  experiments,  were  truncated. 
The  improvement  in  subjective  quality  is  apparent  from  the  photographs. 


9.  Summary 

A  relatively  simple  technique  for  extrapolating  the  values  of  missing 
transform  coefficients  in  transform  signal  processing  systems  has  been 
demonstrated.  The  technique  possesses  the  interesting  feature  that  .t  is 
an  entirely  aposteriori  method.  Thus,  in  a  transform  data  coding  system 
it  is  possible  to  obtain  an  increase  in  data  reconstruction  accuracy  by  pro¬ 
cessing  at  the  decoder,  without  affecting  the  structure  of  the  coder.  Thus, 
in  a  data  coding  system,  it  may  be  worthwhile  to  utilize  a  transform  with 
a  relatively  poor  energy  compaction  property,  but  one  possessing  an 
efficient  computational  algorithm  or  implementation  technology,  and  still 

achieve  a  desired  level  of  system  performance  by  spectrum  extrapolation 
at  the  decoder. 

Several  interesting  extensions  of  the  spectrum  extrapolation  proce 
remain  to  be  explored.  First,  in  many  application,  there  is  a  priori  know¬ 
ledge  available  as  to  the  structure  of  the  reconstructed  data  which  should  be 
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used  in  the  estimation  of  transform  components.  For  example,  in  image 
coding  it  is  known  that  the  image  reconstruction  cannot  possess  negative 
points  or  points  of  extremely  large  magnitude.  A  second  extension  in¬ 
volves  discovery  of  a  method  for  interpolating  the  values  n{  quantized 
transform  coefficients  for  digital  transform  data  coding  systems  in  which 
gross  quantization  is  applied  to  achieve  a  bandwidth  compression.  Finally, 
the  development  of  recursive  estimation  algorithms  to  be  used  in  conjunction 
with  pipeline  transform  generators  would  be  of  interest. 
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This  note  describes  additions  and  modifications  to  a  signal  processing 
interpreter  which  was  developed  as  part  of  the  ARPA  sponsored  program  for  image 
transmission  via  spread  spectrum  links.  The  signal  processing  interpreter 
SPIN3  is  an  interactive  program  for  use  at  a  time-sharing  demand  terminal.  It 
provides  the  user  with  the  equivalent  of  a  calculator  designed  to  perform 
signal  processing  operations  and  provides  the  software  equivalent  of  a  large 
number  of  modules  for  breadboarding  a  complete  signal  processing  system. 

This  interpreter  is  meant  to  be  used  by  engineers  and  scientists  who  are 
familiar  with  signal  processing,  but  who  may  have  no  knowledge  of  programming. 

The  signal  processing  interpreter  is  particularly  useful  for  the  rapid 
investigation  of  systems  whose  complexity  precludes  a  complete  analytic  study, 
and  whose  utilization  of  new  components  may  make  hardware  breadboarding 
undesirable  because  of  cost  and  procurement  time  limitations. 

The  present  note  describes  software  additions  and  modifications  generated 
during  the  period  July  through  September  1973  and  is  meant  to  be  used  in 
conjunction  with  the  Signal  Processing  Interpreter  Preliminary  Description 
and  User's  Guide  [1  ]. 
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NEK  corners  AND  ICMLY  IfflEfNIED  COMWIDS 

quantizer 

corrfct 

noise  generator 

add  . 

clear 

clear  x 

roll  up 

roll  down 

print  all  commands 

scale 

reciprocal 

real  part 

imag  part 

expand 

decimate 

right  circular  shift 


MODIFICATIONS 


1.  A  re-read  has  been  provided  when  read  errors  occur  while  entering  commands 
or  inputting  a  vector.  The  effect  of  this  is  that  noise  on  the  telephone 
link  or  user  format  errors  will  now  only  result  in  an  error  message  and 

a  request  to  re-enter  the  command  or  data  point. 

2.  Previous  input/output  inconsistencies  due  to  differences  between  Fortran 
and  Algol  contentions  have  been  eliminated. 

a.  All  vector  indices  are  now  consistent;  they  run  from  0  to  N-l. 
These  indices  will  appear  in  listings  of  vectors,  plotting, 
and  inputting. 

b.  The  implementation  of  the  command  CONSTANT  FUNCTION  has  been 
rewritten  to  use  a  comma  as  the  delimiter  between  the  real  and 
imaginary  parts.  This  convention  is  now  standard  for  all 
complex  inputs. 

3.  The  output  printing  has  been  cleaned  up  to  minimize  unnecessary  line  feeds 
and  printing  of  redundant  information. 


DETAILED  DESCRIPTION  OF  W  NEW  COIWNDS  AND  NEWLY  IFPlOCMIED  OTH NDS 


QUANTIZER  -  requests  the  user  to  enter  an  integer.  The  x  vector  is  scaled 
so  that  the  largest  absolute  value  of  the  real  or  imaginary  part  of  x  is  set 
equal  to  the  specified  integer.  The  real  and  imaginary  parts  are  then  rounded 
to  the  nearest  integer. 

CORRECT  -  permits  the  user  to  change  the  value  of  a  single  component  of  a 
vector.  The  user  is  asked  to  specify  the  index  of  the  point  to  be  changed, 
and  is  then  asked  to  specify  the  new  value  in  the  format  real  part,  imaginary 
part. 

NOISE  GENERATOR  -  generates  a  zero-mean  random  vector  using  a  multiplicative 
congruential  method  described  by  Kruska  [2].  The  user  supplies  a  value  for  the 
desired  variance  and  an  integer  in  the  range  0-8191  for  the  desired  initializa¬ 
tion  point. 

ADD  ■  pads  the  shorter  of  the  X  and  Y  vectors  with  zeros  to  the  length  of  the 
longer  one.  The  X  vector  is  replaced  by  the  sum  of  X  and  Y,  and  the  stack  is 
lowered. 

"  sets  vectors  to  zero,  with  all  dimensions  set  equal  to  1. 

n  far  y  -  sets  the  X_vector  equal  to  zero,  and  sets  the  X  dimension  qual 
to  one. 

ROU  IJP  -  performs  a  cyclic  permutation  of  the  vectors  X,  Y,  Z,  T  and  a 

similar  permutation  of  their  dimensions  NX,  NY,  NZ,  NT.  Y  is  replaced  by  the 
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old  value  of  X,  etc. 


RQI  l  DOWN  -  performs  the  inverse  of  rqu  HP- 

°R1nt  all  comms  -  causes  the  computer  to  print  back  each  user-entered 
command  which  is  recognized.  This  command  is  used  primarily  in  running  batch 
jobs,  to  provide  a  printout  which  is  readable  without  reference  to  the  card 
deck  which  was  used  for  submission  of  the  job. 

SCALE  -  requests  the  user  to  enter  a  complex  number.  The  X  vector  is 
multiplied  by  the  number,  and  the  result  is  placed  in  X. 

RFCIPROCAI  -  replaces  every  point  in  the  X  vector  by  its  reciprocal. 

REAI  PART  -  replaces  the  X  vector  by  its  real  part. 

I MAG  PART  -  replaces  the  X  vector  by  its  imaginary  part,  i.e.,  1.2, 3. 4  is 

replaced  by  3. 4,0.0. 

EXPAND  "  increases  the  length  of  the  X  vector  by  inserting  zeros  between  the 
original  data  points.  The  '^er  is  requested  to  enter  an  integer  for  the 
expansion  factor.  If  the  integer  k  is  entered,  then  k-1  zeros  are  placed 
between  each  pair  of  data  points  and  after  the  last  point. 

DECIMATE  -  decreases  the  length  of  the  X  vector  by  resampling.  The  user  is 
requested  to  enter  a  positive  integer  for  the  decimation  factor.  If  the  integer 
k  is  entered,  then  every  kth  point  is  saved. 

RIGHT  CIRCULAR  SHIFT  -  produces  a  circular  shift  of  the  X  vector.  The  user 
is  asked  to  enter  an  integer  for  the  shift  value.  If  the  integer  is  positive, 
it  causes  a  right  circular  shift.  If  the  integer  is  negative,  it  causes  a  left 
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circular  shift,  since  a  circular  shift  by  a  multiple  of  NX  is  equivalent  to 

a  Shift  of  zero,  a  circular  shift  of  -k  is  equivalent  to  a  circular  shift  of 

NX-k,  so  the  distinction  between  right  circular  shifts  and  left  circular  shifts 
is  purely  conceptual. 
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SIGNAL  TRANSFORMATION  WITH  RECIRCULATING  SAW  DELAY  LINES* 

by 
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Stanford,  California 

ABSTRACT 

A  system  is  described  for  transforming  data  by  switching  it  into  a 
bank  of  recirculating  delay  lines  in  one  time  sequence  and  switching  it 
out  in  another  sequence.  This  is  applied  to  row-to-coluran  conversion 
when  the  input  data  is  in  the  form  of  an  N  x  N  matrix.  An  experimental 
demonstration  is  shown  which  uses  SAW  wrap-around  delay  line  channels, 
transforming  a  3  y  3  matrix  of  rf  pulses  of  16  psec  length  and  52  MHz 
carrier  frequency. 


* 

The  work  reported  in  this  paper  was  supported  by  the  U.S.  Office  o 
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SIGNAL  TRANSFORMATION  WITH  RECIRCULATING  SAW  DELAY  LINES 

by 

C>.  M.  Fortunko  and  h.  J.  Shaw 

I .  INTRODUCTION 

Data  processing  in  real  time  can  be  accomplished  by  storing  incoming 
data  in  banks  of  recirculating  surface  wave  delay  lines  and  then  reading 
out  these  data  in  such  a  way  as  to  transform  them  in  a  desired  fashion. 

We  consider  here  some  of  the  properties  of  such  systems  and  describe  the 
operation  of  a  simple  prototype  system,  operating  at  52  MHz  carrier  fre¬ 
quency,  which  performs  an  operation  equivalent  to  interchanging  the  rows 
and  columns  of  a  matrix. 

The  basic  approach  to  the  use  of  delay  lines  for  such  purposes  is 
illustrated  in  Fig.  1.  The  essential  item  is  a  bank  of  identical 
recirculating  surface  wave  delay  line  loops.  Input  data  fed  in  at  the  input 
terminal  of  the  system  are  distributed  among  the  delay  lines  in  some  desired 
fashion  by  means  of  the  WRITE  circuitry.  The  data  enter  the  delay  lines 
through  interdigital  transducers  indicated  as  arrows  in  the  figure.  The  READ 
circuitry  extracts  the  data  from  the  delay  lines  by  switching  from  one 
delry  line  to  another  in  a  time  sequence  which  is  different  from  that  used 
for  writing,  to  accomplish  the  required  transformation  of  the  data  train. 

The  row-to-column  converter  can  be  used  to  perform  the  intermediate 
serial- to -parallel  conversion  which  is  required  to  implement  decomposable 
two-dimensional  transforms  (e.g.  Fourier,  Hadamard).  It  can  also  be  used 
to  interface  two  system  with  different  speeds. 


II.  ROW-TO-COLUMN  CONVERSION  USING 
CLOSED-LOOP  SURFACE  ACOUSTIC  WAVE  DELAY  LINES 


He  consider  now  an  example  of  data  transformation.  Let  the  input 
data  stream  be  An  •  •  •  *1N  .  .  .  AJ(1  .  .  .  A^,  which  is  ordered 
according  to  rows  of  the  N  X  N  matrix  in  Fig.  2.  He  now  show  that  it  is 
posaibie  to  devise  WRITE  and  READ  aigorlthms  which  will  convert  this 
into  a  sequence  ordered  by  columns  of  the  matrix,  in  real  time,  using  only 
a  single  transducer  on  each  delay  line. 

The  transducers  of  Fig.  1  can  be  addressed  at  either  of  two  rates, 
once  per  row  or  once  per  bit.  To  obtain  the  row-to-column  conversion  one 
of  these  rates  is  applied  to  the  WRITE  function  and  the  other  to  the  READ 
function.  Either  rate  may  be  applied  to  either  function  giving  two  algo¬ 
rithms  for  the  row-to-column  operation.  Figures  3(a)  and  (b)  illustrate 
these  two  algorithms.  In  these  diagrams  successive  horizontal  lines 
show  the  signals  seen  by  the  transducer  in  successive  delay  lines  as  a 
function  of  time  reading  to  the  right.  The  squares  represent  the  signals 
introduced  into  the  transducers  during  the  WRITE  operation,  and  the 
quantities  between  curved  brackets  are  the  signals  extracted  from  the 

transducers  in  the  READ  operation.  The  dashes  indicate  time  slots 
containing  no  signal. 


Figure  3(a)  is  for  the  case  in  shich  the  WRITE  function  is  performed 
at  the  row  rate  and  the  READ  function  at  the  bit  rate.  Each  row  of  Fig.  2 
is  fed  intact  into  a  separate  delay  line.  As  indicated  by  the  diagonal 
arrays  of  curved  brackets,  the  sequence  of  gating  on  the  transducers  for 


readout  is  a  uniform,  periodic  scan,  1,  2  ...  N  ,  1,2 
sequentially  from  line  1  to  line  N,  in  repetitive  cycles 


The  final  output 
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signal  is  shown  at  the  bottom  of  Fig.  3v.i).  This  output  sequence  is  a  real¬ 
time  sequence  free  of  gaps,  ordered  according  to  columns  (time  reversed) 
of  the  matrix  of  Fig.  2. 

Figure  3(b)  represents  the  alternate  algorithm  in  which  the  WRITE 
and  READ  switching  rates  have  been  interchanged.  For  any  row  of  input 
data  of  Fig.  2,  the  first  bit  enters  the  first  storage  loop,  the  second 
bit  enters  the  second  loop,  and  so  forth  until  the  entire  row  has  been 
stored,  one  bit  in  each  delay  line.  During  the  READ  cycle  the  entire 
-'ontents  of  each  delay  line  are  sampled,  one  row  at  a  time.  The  sequence 
(column  output)  shown  at  the  bottom  of  Fig.  3(b)  is  seen  to  be  the  same 
as  for  Fig.  3(a),  ordered  according  to  the  columns  (time  reversed'*  of  the 
matrix  of  Fig.  2.  We  see  that  the  algorithms  of  Figs,  rva1*  and  (b' 
result  in  output  sequences  whose  rows  are  columns  of  the  input  matrix 
as  desired.  However,  both  the  rows  and  the  elements  within  these  rows 
appear  in  the  opposite  time  sequence  to  the  corresponding  columns  of  Fig.  2. 
That  is,  the  output  sequences  in  Figs.  3(a)  and  (b)  are  time  reversed.  For 
cases  of  direct  concern  here,  this  is  not  of  consequence.  If  it  is  not 
desired,  however,  it  can  be  corrected  by  time-reversing  the  input  signal 
by  rows  before  feeding  into  the  delay  line  loops.  Figure  3(c)  shows  this 
behavior.  The  signals  within  each  row  have  been  time-reversed  before 
insertion  into  the  delay  line  loops.  This  time  reversal  could  be  done 
in  a  device  having  a  time  delay  of  N  clock  intervals,  rather  than  N  X  N 
as  would  be  required  if  row-to-row  time  reversal  were  also  involved.  A 
parametric  acoustic  time  reversal  could  be  employed  for  this  operation. 
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Fig.  2- -Matrix  Notation. 
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* -Gating  sequences  for  row-to-column  conversion. 


